1 SETUP

1.1 Open libraries

library(gtsummary)
library(dplyr)
library(tidyverse)
library(survival)
library(survminer)
library(lubridate)
library(readxl)
library(janitor)
library(ggplot2)
library(tidycmprsk)
library(ggsurvfit)
library(tidymodels) 
library(sjPlot)
library(patchwork)
library(bstfun)
library(ComplexUpset)

library(writexl)
library(rms)

library(timeROC)
library(rms)
library(pec)
library(prodlim)
library(riskRegression)

library(plotly)
library(ggtext)

library(WeightIt)
library(cobalt)
library(survey)

2 DATA WRANGLING

2.1 Read in MM datasets

db_Kansas = read_excel(path = "/Users/aportugu/Desktop/Elranatamab/Analysis/Data/bsAb Consortium Data Template 2.11.2025 - due 5-1-25 (de-identified).xlsx", sheet = "bsAb_DATA")%>%
  filter(`bsAb name` %in% c("1","3"))%>%
  clean_names()%>%
  filter(!is.na(age_at_bs_ab_initiation)) %>%
  mutate(
    patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate = as.character(patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate),
    baseline_anc = as.double(baseline_anc),
    baseline_alc = as.double(baseline_alc),
    alc_day_8 = as.double(alc_day_8),
    alc_day_15 = as.double(alc_day_15),
    alc_day_22 = as.double(alc_day_22),
    alc_c2d1 = as.double(alc_c2d1),
    anc_nadir_value_during_first_6_cycles = as.double(anc_nadir_value_during_first_6_cycles),
    cr_cl_prior_to_step_up_dose_number_1 = as.double(cr_cl_prior_to_step_up_dose_number_1)
  )

# Write to Excel
#write_xlsx(db_Kansas, "Kansas_elra_data.xlsx")

db_CC = read_excel(path = "/Users/aportugu/Desktop/Elranatamab/Analysis/Data/CC bsAb Consortium Submission CCF 6-11-25.xlsx", sheet = "bsAb_DATA")%>%
  filter(`bsAb name` %in% c("1","3"))%>%
  clean_names()%>%
  mutate(
    patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate = as.character(patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate),
    baseline_anc = as.double(baseline_anc),
    baseline_alc = as.double(baseline_alc),
    alc_day_8 = as.double(alc_day_8),
    alc_day_15 = as.double(alc_day_15),
    alc_day_22 = as.double(alc_day_22),
    alc_c2d1 = as.double(alc_c2d1),
    anc_nadir_value_during_first_6_cycles = as.double(anc_nadir_value_during_first_6_cycles),
    cr_cl_prior_to_step_up_dose_number_1 = as.double(cr_cl_prior_to_step_up_dose_number_1)
  )

# Write to Excel
#write_xlsx(db_CC, "CC_elra_data.xlsx")

db_HUMC = read_excel(path = "/Users/aportugu/Desktop/Elranatamab/Analysis/Data/BsAb Consortium Data HUMC.xlsx", sheet = "bsAb_DATA")%>%
  filter(`bsAb name` %in% c("1","3"))%>%
  clean_names()%>%
  mutate(
    patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate = as.character(patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate),
    baseline_anc = as.double(baseline_anc),
    baseline_alc = as.double(baseline_alc),
    alc_day_8 = as.double(alc_day_8),
    alc_day_15 = as.double(alc_day_15),
    alc_day_22 = as.double(alc_day_22),
    alc_c2d1 = as.double(alc_c2d1),
    anc_nadir_value_during_first_6_cycles = as.double(anc_nadir_value_during_first_6_cycles),
    cr_cl_prior_to_step_up_dose_number_1 = as.double(cr_cl_prior_to_step_up_dose_number_1)
  )

# Write to Excel
#write_xlsx(db_HUMC, "HUMC_elra_data.xlsx")

db_MCC = read_excel(path = "/Users/aportugu/Desktop/Elranatamab/Analysis/Data/bsAb Consortium Data MCC-deidentified.xlsx", sheet = "bsAb_DATA")%>%
  filter(`bsAb name` %in% c("1","3"))%>%
  clean_names()%>%
  mutate(
    patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate = as.character(patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate),
    anc_nadir_value_during_first_6_cycles = as.double(anc_nadir_value_during_first_6_cycles),
    baseline_anc = as.double(baseline_anc),
    baseline_alc = as.double(baseline_alc),
    cr_cl_prior_to_step_up_dose_number_1 = as.double(cr_cl_prior_to_step_up_dose_number_1)
  )


db_Huntsman = read_excel(path = "/Users/aportugu/Desktop/Elranatamab/Analysis/Data/Utah_Huntsman_bsAb Consortium Data FINAL_5.4.25.xlsx", sheet = "bsAb_DATA")%>%
  filter(`bsAb name` %in% c("1","3"))%>%
  clean_names()%>%
  mutate(
    patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate = as.character(patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate),
    baseline_anc = as.double(baseline_anc)/1000,
    baseline_alc = as.double(baseline_alc)/1000,
    anc_nadir_value_during_first_6_cycles = as.double(anc_nadir_value_during_first_6_cycles),
    cr_cl_prior_to_step_up_dose_number_1 = as.double(cr_cl_prior_to_step_up_dose_number_1)
  )

# Write to Excel
#write_xlsx(db_Huntsman, "Huntsman_elra_data.xlsx")


db_FHCC = read_excel(path = "/Users/aportugu/Desktop/Elranatamab/Analysis/Data/FHCC Bispecifics (5.5.2025).xlsx", sheet = "bsAb_DATA")%>%
  filter(`bsAb name` %in% c("1","3"))%>%
  clean_names()%>%
  mutate(
    patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate = as.character(patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate),
    date_of_infection_80 = ymd(date_of_infection_80),
    baseline_anc = as.double(baseline_anc),
    cr_cl_prior_to_step_up_dose_number_1 = as.double(cr_cl_prior_to_step_up_dose_number_1)
  )


db_UTSW = read_excel(path = "/Users/aportugu/Desktop/Elranatamab/Analysis/Data/bsAb Consortium Data Template 2.11.2025 1.xlsx", sheet = "bsAb_DATA")%>%
  filter(`bsAb name` %in% c("1","3")) %>%
  clean_names()%>%
  mutate(
    patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate = as.character(patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate),
    baseline_anc = as.double(baseline_anc),
    baseline_alc = as.double(baseline_alc),
    anc_nadir_value_during_first_6_cycles = as.double(anc_nadir_value_during_first_6_cycles)
  )

# Write to Excel
#write_xlsx(db_UTSW, "UTSW_elra_data.xlsx")


db_MDACC = read_excel(path = "/Users/aportugu/Desktop/Elranatamab/Analysis/Data/MDACC BITE ASH 2025.xlsx", sheet = "bsAb_DATA")%>%
  filter(`bsAb name` %in% c("1","3"))%>%
  clean_names()%>%
  mutate(
    patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate = as.character(patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate),
    baseline_anc = as.double(baseline_anc),
    baseline_alc = as.double(baseline_alc),
    anc_nadir_value_during_first_6_cycles = as.double(anc_nadir_value_during_first_6_cycles),
    cr_cl_prior_to_step_up_dose_number_1 = as.double(cr_cl_prior_to_step_up_dose_number_1),
    total_number_paraskelatal_plasmacytomas = as.double(total_number_paraskelatal_plasmacytomas)
  )

# Write to Excel
#write_xlsx(db_MDACC, "MDACC_elra_data.xlsx")


db_Roswell = read_excel(path = "/Users/aportugu/Desktop/Elranatamab/Analysis/Data/bsAb Consortium Data Template 5.8.2025_NN.xlsx", sheet = "bsAb_DATA")%>%
  filter(`bsAb name` %in% c("1","3"))%>%
  clean_names()%>%
  mutate(
    patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate = as.character(patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate),
    baseline_anc = as.double(baseline_anc),
    baseline_alc = as.double(baseline_alc),
    anc_nadir_value_during_first_6_cycles = as.double(anc_nadir_value_during_first_6_cycles),
    cr_cl_prior_to_step_up_dose_number_1 = as.double(cr_cl_prior_to_step_up_dose_number_1),
    total_number_paraskelatal_plasmacytomas = as.double(total_number_paraskelatal_plasmacytomas)
  )

# Write to Excel
#write_xlsx(db_Roswell, "Roswell_elra_data.xlsx")

# Combine all into a single dataframe
combined_db <- bind_rows(
  list(
    HUMC = db_HUMC,
    MCC = db_MCC,
    Huntsman = db_Huntsman,
    FHCC = db_FHCC,
    UTSW = db_UTSW,
    MDACC = db_MDACC,
    Roswell = db_Roswell,
    CC = db_CC,
    Kansas = db_Kansas
  ),
  .id = "site"
) %>%
  mutate(
    prior_bcma_directed_therapy_yes_no = ifelse(!is.na(most_recent_bcma_agent_name),1,0),
    bs_ab_name = case_when(
      bs_ab_name == 1 ~ "tec",
      bs_ab_name == 3 ~ "elra",
      .default = NULL
    ),
    max_crs_grade_0_5 = ifelse(is.na(max_crs_grade_0_5), 0,max_crs_grade_0_5), ## Set NAs to 0
    max_icans_grade = ifelse(is.na(max_icans_grade), 0, max_icans_grade ) ## Set NAs to 0
  )

# Convert all datetime columns to Date
combined_db_clean <- combined_db %>%
  mutate(across(where(lubridate::is.POSIXt), as.Date))

# Write to Excel
write_xlsx(combined_db_clean, "elra_tec_data.xlsx")


# Write to Excel
# write_xlsx(combined_db_clean %>% filter(prior_bcma_directed_therapy_yes_no != prior_bcma_directed_therapy), "elra_data_discordant.xlsx")

## Several institutions do not have any cases of elranatamab

# db_UABMC = read_excel(path = "Data/bsAb Consortium Data 5.8.xlsx", sheet = "bsAb_DATA")%>%
#   filter(`bsAb name` == "3")%>%
#   clean_names()

# db_Stanford = read_excel(path = "Data/De-identified Stanford bsAb Consortium Data Template 4.30.2025.xlsx", sheet = "bsAb_DATA")%>%
#   filter(`bsAb name` == "3")%>%
#   clean_names()

# db_LCI = read_excel(path = "Data/bsAb Consortium Data LCI.xlsx", sheet = "bsAb_DATA")%>%
#   filter(`bsAb name` == "3")%>%
#   clean_names()

# db_COH = read_excel(path = "Data/bsAb Consortium Data Template 2.11.2025_COHv2.xlsx", sheet = "bsAb_DATA")%>%
#   filter(`bsAb name` == "3")%>%
#   clean_names()

# db_MUSC = read_excel(path = "Data/USMM bsAb Consortium Data Template 2.11.2025.xlsx", sheet = "bsAb_DATA")%>%
#   filter(`bsAb name` == "3")%>%
#   clean_names()

# db_Mayo = read_excel(path = "Data/Bispecific_consortium_20250521 NoPHI.xlsx", sheet = "bsAb_DATA") %>%
#   filter(`bsAb name` == 3)%>%
#   clean_names()

2.2 Calculate CAR-HEMATOTOX score and risk level

combined_db <- combined_db %>% 
  mutate(
    # baseline_anc = as.numeric(baseline_anc),
    # baseline_hgb,
    # baseline_pl_ts,
    # baseline_crp_prior_to_bs_ab_initiation,
    # baseline_ferritin_prior_to_bs_ab_initiation,
    HT_Score = 
      (baseline_anc <= 1.2) +
      (baseline_hgb <= 9.0) +
      (baseline_pl_ts >= 76 & baseline_pl_ts <= 175) +
      (baseline_crp_prior_to_bs_ab_initiation >= 3.0) +
      (baseline_ferritin_prior_to_bs_ab_initiation >= 650 & baseline_ferritin_prior_to_bs_ab_initiation <= 2000) +
      (baseline_pl_ts <= 75) + ## Add an extra point for Plt <=75
      (baseline_ferritin_prior_to_bs_ab_initiation > 2000), ## Add an extra point for ferritin >2000
    
    # Classify the risk based on the score
    HT_Risk = factor(
      case_when(
        HT_Score >= 2 ~ "HThigh",
        baseline_pl_ts <=75 ~ "HThigh", ## Equivalent to at least 2 points
        baseline_ferritin_prior_to_bs_ab_initiation > 2000 ~ "HThigh", ## Equivalent to at least 2 points
        (baseline_anc <= 1.2) + (baseline_hgb <= 9.0) + (baseline_pl_ts >= 76 & baseline_pl_ts <= 175) >=2 ~ "HThigh",
        HT_Score <2 ~ "HTlow",
        .default = "NA"
      ),
      levels = c("HTlow", "HThigh")
    )
)
write_xlsx(combined_db %>% select(site, patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate, baseline_anc, baseline_hgb, baseline_pl_ts, baseline_crp_prior_to_bs_ab_initiation, baseline_ferritin_prior_to_bs_ab_initiation, HT_Score, HT_Risk), "HT_score.xlsx")

2.3 Identify severe infections

combined_db <- combined_db %>% 
  mutate(
    severe_infection = ifelse(infection_severity_73 %in% 1 | infection_severity_76 %in% 1 | infection_severity_79 %in% 1| infection_severity_82 %in% 1,1,0),
    first_severe_infection_date = case_when(
      infection_severity_73 %in% 1 ~ date_of_infection_74,
      infection_severity_76 %in% 1 ~ date_of_infection_77,
      infection_severity_79 %in% 1 ~ date_of_infection_80,
      infection_severity_82 %in% 1 ~ date_of_infection_83,
      .default = NULL
    ),
    first_severe_infection_type = case_when(
      infection_severity_73 %in% 1 ~ number_infection_1_bacterial_viral_fungal,
      infection_severity_76 %in% 1 ~ number_infection_2_bacterial_viral_fungal,
      infection_severity_79 %in% 1 ~ number_infection_3_bacterial_viral_fungal,
      infection_severity_82 %in% 1 ~ number_infection_4_bacterial_viral_fungal,
      .default = NULL
    )
    
    
  )

3 PROPENSITY SCORES

3.1 Regression analysis for treatment allocation

theme_gtsummary_compact()
## Setting theme "Compact"
variables <- c("elra", "Age", "Female sex", "ECOG 2+", "EMD", "Hgb", "LDH", "Prior BCMA")

data1 <- combined_db %>% 
  drop_na(age_at_bs_ab_initiation, ecog_ps_at_bs_ab_initiation_0_5, just_prior_to_bs_ab_initiation_is_there_true_emd_if_none_move_to_column_gk, baseline_hgb, baseline_ldh_prior_to_bs_ab_initiation) %>%
  transmute(
    Age = age_at_bs_ab_initiation,
    "Female sex" = sex_0_male_1_female,
    elra = ifelse(bs_ab_name == "elra", 1, 0),
    "ECOG 2+" = ifelse(ecog_ps_at_bs_ab_initiation_0_5 >=2,1,0),
    EMD = ifelse(just_prior_to_bs_ab_initiation_is_there_true_emd_if_none_move_to_column_gk == 2,1,0),
    Hgb = as.double(baseline_hgb),
    LDH = as.double(baseline_ldh_prior_to_bs_ab_initiation),
    "Prior BCMA" = prior_bcma_directed_therapy_yes_no
  )

uv_tab <- tbl_uvregression(
    data1[c(variables)],
    method = glm,
    y = "elra",
    method.args = list(family = binomial),
    exponentiate = TRUE
  ) %>%
  sort_p()


mv_tab<-glm(elra ~ Age + `Female sex` + `ECOG 2+` + EMD + Hgb + LDH + `Prior BCMA`, data=data1, family=binomial) %>%
  tbl_regression(exponentiate=TRUE)

tbl_merge(list(uv_tab, mv_tab), tab_spanner = c("**Univariate**", "**Multivariable**"))
Characteristic
Univariate
Multivariable
N OR1 95% CI1 p-value OR1 95% CI1 p-value
Hgb 371 1.12 1.01, 1.25 0.029 1.14 1.02, 1.28 0.026
EMD 371 1.69 0.95, 2.96 0.069 1.71 0.96, 3.04 0.068
Age 371 1.02 1.00, 1.04 0.087 1.01 0.99, 1.04 0.3
Female sex 371 1.43 0.93, 2.21 0.11 1.43 0.92, 2.24 0.11
Prior BCMA 371 0.76 0.49, 1.17 0.2 0.78 0.50, 1.23 0.3
ECOG 2+ 371 1.19 0.75, 1.88 0.4 1.27 0.78, 2.07 0.3
LDH 371 1.00 1.00, 1.00 0.8 1.00 1.00, 1.00 >0.9
1 OR = Odds Ratio, CI = Confidence Interval

3.2 IPTW

IPTW.data <- combined_db %>% 
  drop_na(age_at_bs_ab_initiation, ecog_ps_at_bs_ab_initiation_0_5, just_prior_to_bs_ab_initiation_is_there_true_emd_if_none_move_to_column_gk, baseline_hgb, baseline_ldh_prior_to_bs_ab_initiation) %>%
  mutate(
    Age = age_at_bs_ab_initiation,
    "Female sex" = sex_0_male_1_female,
    bsAb = bs_ab_name,
    "Received full dose" = ifelse(!is.na(date_of_first_bs_ab_full_dose),1,0),
    "ECOG" = case_when(
      ecog_ps_at_bs_ab_initiation_0_5 %in% c("0","1") ~"0-1",
      ecog_ps_at_bs_ab_initiation_0_5 ==2 ~"2",
      ecog_ps_at_bs_ab_initiation_0_5 >=3 ~"3+",
      .default = NULL),
    EMD = ifelse(just_prior_to_bs_ab_initiation_is_there_true_emd_if_none_move_to_column_gk == 2,1,0),
    Hgb = as.double(baseline_hgb),
    LDH = as.double(baseline_ldh_prior_to_bs_ab_initiation),
    "Prior BCMA" = prior_bcma_directed_therapy_yes_no
  )

## Use the WeightIt function to generate propensity scores and weights
W.out <- weightit(bsAb ~ Age + `ECOG` + EMD + Hgb + LDH + `Prior BCMA`,
                  data=IPTW.data, 
                  method = "glm",    ## Methods: glm, gbm, cbps, npcbps, ebal, optweight, super, bart, energy
                  estimand = "ATE",
                  stabilize = TRUE)

bal.tab(W.out, stats = c("m", "v"), thresholds = c(m=0.25), binary = "std")
## Balance Measures
##                Type Diff.Adj     M.Threshold V.Ratio.Adj
## prop.score Distance   0.0050 Balanced, <0.25      0.9733
## Age         Contin.   0.0090 Balanced, <0.25      0.9595
## ECOG_0-1     Binary  -0.0121 Balanced, <0.25           .
## ECOG_2       Binary   0.0152 Balanced, <0.25           .
## ECOG_3+      Binary  -0.0027 Balanced, <0.25           .
## EMD          Binary  -0.0034 Balanced, <0.25           .
## Hgb         Contin.  -0.0186 Balanced, <0.25      1.0974
## LDH         Contin.  -0.0276 Balanced, <0.25      1.0961
## Prior BCMA   Binary  -0.0119 Balanced, <0.25           .
## 
## Balance tally for mean differences
##                     count
## Balanced, <0.25         9
## Not Balanced, >0.25     0
## 
## Variable with the greatest mean difference
##  Variable Diff.Adj     M.Threshold
##       LDH  -0.0276 Balanced, <0.25
## 
## Effective sample sizes
##              elra    tec
## Unadjusted 123.   248.  
## Adjusted   113.38 242.39
summary(W.out)
##                   Summary of weights
## 
## - Weight ranges:
## 
##         Min                                  Max
## elra 0.4745 |---------------------------| 2.3134
## tec  0.7811      |-------------|          1.7418
## 
## - Units with the 5 most extreme weights by group:
##                                         
##          123    102     81     91      9
##  elra 1.5458 1.5896  1.653 1.7607 2.3134
##          169     71    274    258    271
##   tec 1.4257 1.4453 1.5215 1.6783 1.7418
## 
## - Weight statistics:
## 
##      Coef of Var   MAD Entropy # Zeros
## elra       0.292 0.222   0.040       0
## tec        0.152 0.115   0.011       0
## 
## - Effective Sample Sizes:
## 
##              elra    tec
## Unweighted 123.   248.  
## Weighted   113.38 242.39
love_plot <- love.plot(W.out,
          binary = "std",
          line = TRUE, 
          abs = TRUE,
          thresholds = c(m=0.2),
          sample.names = c("Unweighted", "IPTW"),
          title = "A)",
          var.order = "unadjusted",
            var.names = c(
              prop.score = "Distance"
          ))+ theme(plot.title = element_text(hjust = 0))
love_plot

weight_df = data.frame(W.out$weights, W.out$treat)

IPTW.data$weights <- W.out$weights
  
IPTW.data$ps <- W.out$ps


write.csv(data1, "IPTW.dataset.csv")

3.3 Scatterplot of weights

Weights_plot <- ggplot(
  IPTW.data, aes(x=ps, y = weights, color = bsAb)) + 
  geom_point(shape = 16) +
  theme_classic() + 
  scale_y_continuous(breaks = seq(0,4 , by = 0.5)) +
  labs(x = "Propensity score", y = "Weights") + 
  ggtitle("B)")
Weights_plot

4 TABLES

4.1 Site contributions

theme_gtsummary_compact()

combined_db %>%
  transmute(
    "Site" = site,
    bs_ab_name,
  ) %>%
  tbl_summary(
    by = bs_ab_name,
    missing = "no",
    sort = list(all_categorical() ~ "frequency"),
    statistic = list(
      all_continuous() ~ "{median} ({p25}, {p75})",
      all_categorical() ~ "{n} ({p}%)"
    )
  ) %>%
  bold_labels()
Characteristic elra
N = 130
1
tec
N = 277
1
Site

    MCC 53 (41%) 72 (26%)
    MDACC 12 (9.2%) 57 (21%)
    CC 20 (15%) 47 (17%)
    UTSW 10 (7.7%) 40 (14%)
    FHCC 7 (5.4%) 24 (8.7%)
    Huntsman 9 (6.9%) 22 (7.9%)
    Kansas 4 (3.1%) 13 (4.7%)
    HUMC 10 (7.7%) 1 (0.4%)
    Roswell 5 (3.8%) 1 (0.4%)
1 n (%)

4.2 Baseline characteristics

theme_gtsummary_compact()

combined_db %>%
  transmute(
    #"Site" = site,
    "Age (years)" = age_at_bs_ab_initiation,
    "Female sex" = sex_0_male_1_female == 1,
    "Race" = factor(case_when(
      is.na(race_white_0_black_1_aapi_2_american_indian_3_other_4_unk_5) ~ NA,
      race_white_0_black_1_aapi_2_american_indian_3_other_4_unk_5 == 0 ~ "White",
      race_white_0_black_1_aapi_2_american_indian_3_other_4_unk_5 == 1 ~ "Black",
      race_white_0_black_1_aapi_2_american_indian_3_other_4_unk_5 == 2 ~ "AAPI",
      race_white_0_black_1_aapi_2_american_indian_3_other_4_unk_5 == 3 ~ "American Indian",
      race_white_0_black_1_aapi_2_american_indian_3_other_4_unk_5 == 4 ~ "Other",
      race_white_0_black_1_aapi_2_american_indian_3_other_4_unk_5 == 5 ~ NA
    ), levels = c("White","Black","AAPI","American Indian","Other",NA)),
    #ecog_ps_at_bs_ab_initiation_0_5,
    "ECOG" = case_when(
      ecog_ps_at_bs_ab_initiation_0_5 == 0 ~ "0",
      ecog_ps_at_bs_ab_initiation_0_5 == 1 ~ "1",
      ecog_ps_at_bs_ab_initiation_0_5 == 2 ~ "2",
      ecog_ps_at_bs_ab_initiation_0_5 >= 3 ~ "3+"
    ),
    "LDH (U/L)" = as.double(baseline_ldh_prior_to_bs_ab_initiation),
    "Hgb (g/dL)" = baseline_hgb,
    "Heavy chain" = factor(case_when(
      myeloma_involved_heavy_chain_ig_g_ig_a_ig_m_ig_d_ig_e_or_blank == "IgG" ~ "IgG",
      myeloma_involved_heavy_chain_ig_g_ig_a_ig_m_ig_d_ig_e_or_blank == "IgM" ~ "IgM",
      myeloma_involved_heavy_chain_ig_g_ig_a_ig_m_ig_d_ig_e_or_blank == "IgA" ~ "IgA",
      myeloma_involved_heavy_chain_ig_g_ig_a_ig_m_ig_d_ig_e_or_blank == "IgD" ~ "IgD",
      .default = "None"
      ), levels = c("IgG", "IgA","IgM","IgD","None")),
    "True EMD" = just_prior_to_bs_ab_initiation_is_there_true_emd_if_none_move_to_column_gk == 2,
    "Oligo or non-secretory" = oligo_or_non_secretory_at_bs_ab_initiation_0_no_1_yes,
    "≥1 HRCA" = ifelse(deletion_17p_prior_to_bs_ab_0_no_1_yes==1 | t_4_14_prior_to_bs_ab_0_no_1_yes == 1 | t_14_16_prior_to_bs_ab_0_no_1_yes == 1 | amp_1_q_only_0_no_1_yes == 1 | gain_amp1q21_prior_to_bs_ab_0_no_1_yes == 1,1,0),
    "del(17p)" = deletion_17p_prior_to_bs_ab_0_no_1_yes,
    "t(4;14)" = t_4_14_prior_to_bs_ab_0_no_1_yes,
    "t(14;16)" = t_14_16_prior_to_bs_ab_0_no_1_yes,
    "gain/amp(1q)" = gain_amp1q21_prior_to_bs_ab_0_no_1_yes,
    "amp(1q)" = amp_1_q_only_0_no_1_yes,
    "Prior LOTs" = number_of_prior_lines_of_therapy_before_bs_ab,
    "Prior ASCT" = prior_auto_sct_no_0_yes_1,
    "Prior CAR-T" = prior_cart_0_no_1_yes,
    "Bortezomib refractory" = bortezomib_status == 1,
    "Lenalidomide refractory" = lenalidomide_status == 1,
    "Pomalidomide refractory" = pomalidomide_status == 1,
    "Anti-CD38 refractory" = daratumumab_or_isatuximab_status == 1,
    "Triple refractory" = triple_class_refractory_i_mi_d_pi_anti_cd38 == 1,
    "Penta refractory" = penta_refractory_2_p_is_2_i_mi_ds_anti_cd38 == 1,
    "Prior BCMA-directed therapy" = !is.na(most_recent_bcma_agent_name),
    "Most recent BCMA" = factor(case_when(
      most_recent_bcma_agent_name %in% c("Abecma", "ABECMA","Abecma on trial","Ide-cel") ~ "CAR-T",
      most_recent_bcma_agent_name %in% c("Carvykti","carvykti", "Carvikty", "Clita-cel", "Cilta-cel") ~ "CAR-T",
      most_recent_bcma_agent_name %in% c("CC98633 CAR-T","ddBCMA") ~ "CAR-T",
      most_recent_bcma_agent_name %in% c("BLENREP", "blenrep", "Belantamab") ~ "ADC",
      most_recent_bcma_agent_name %in% c("Teclistamab") ~ "TCE"
    ), levels = c("CAR-T","TCE", "ADC")),
    "Trial eligible" = eligible_for_corresponding_bs_ab_trial_on_first_dose_use_trial_elig_sheet_to_help_you == 1,
    "Received full dose" = ifelse(!is.na(date_of_first_bs_ab_full_dose),1,0),
     "Baseline ferritin (ng/mL)" = as.double(baseline_ferritin_prior_to_bs_ab_initiation),
     "Baseline Hgb" = as.double(baseline_hgb),
     "Baseline LDH" = as.double(baseline_ldh_prior_to_bs_ab_initiation),
    "Duration of prior BCMA" = as.numeric(ymd(date_of_last_exposure_to_prior_bcma_agent) - ymd(date_of_first_exposure_to_prior_bcma_agent)),
    bs_ab_name
  ) %>%
  tbl_summary(
    by = bs_ab_name,
    missing = "no",
    #sort = list(all_categorical() ~ "frequency"),
    statistic = list(
      all_continuous() ~ "{median} ({p25}, {p75})",
      all_categorical() ~ "{n} ({p}%)",
      "Age (years)" ~ "{median} ({min} to {max})"
    )
  ) %>%
  add_p() %>%
  bold_labels() %>%
  add_variable_grouping(
    "High-risk cytogenetic abnormality at any time" = c("del(17p)","t(4;14)", "t(14;16)", "gain/amp(1q)", "amp(1q)"),
    "Refractory status" = c("Bortezomib refractory", "Lenalidomide refractory", "Pomalidomide refractory", "Anti-CD38 refractory")
  )
Characteristic elra
N = 130
1
tec
N = 277
1
p-value2
Age (years) 71 (39 to 95) 69 (35 to 88) 0.085
Female sex 74 (57%) 128 (46%) 0.044
Race

0.14
    White 102 (82%) 211 (76%)
    Black 16 (13%) 54 (20%)
    AAPI 3 (2.4%) 7 (2.5%)
    American Indian 2 (1.6%) 0 (0%)
    Other 2 (1.6%) 4 (1.4%)
ECOG

0.5
    0 22 (17%) 52 (20%)
    1 60 (48%) 130 (49%)
    2 30 (24%) 65 (25%)
    3+ 14 (11%) 18 (6.8%)
LDH (U/L) 223 (177, 277) 208 (166, 287) 0.4
Hgb (g/dL) 10.05 (8.70, 11.70) 9.60 (8.30, 11.30) 0.020
Heavy chain

<0.001
    IgG 75 (58%) 122 (44%)
    IgA 28 (22%) 47 (17%)
    IgM 2 (1.5%) 1 (0.4%)
    IgD 2 (1.5%) 1 (0.4%)
    None 23 (18%) 106 (38%)
True EMD 28 (22%) 36 (13%) 0.028
Oligo or non-secretory

0.3
    0 116 (91%) 259 (94%)
    1 12 (9.4%) 15 (5.4%)
    2 0 (0%) 2 (0.7%)
≥1 HRCA 85 (69%) 164 (63%) 0.2
High-risk cytogenetic abnormality at any time


    del(17p) 31 (25%) 70 (27%) 0.6
    t(4;14) 17 (13%) 34 (13%) >0.9
    t(14;16) 8 (6.3%) 15 (5.8%) 0.8
    gain/amp(1q) 70 (56%) 128 (49%) 0.2
    amp(1q) 9 (7.9%) 18 (6.9%) 0.7
Prior LOTs 6.00 (4.00, 8.00) 6.00 (4.00, 8.00) 0.8
Prior ASCT 66 (51%) 163 (59%) 0.13
Prior CAR-T 54 (42%) 103 (37%) 0.4
Refractory status


    Bortezomib refractory 87 (67%) 226 (82%) <0.001
    Lenalidomide refractory 110 (85%) 238 (86%) 0.7
    Pomalidomide refractory 96 (74%) 226 (82%) 0.062
    Anti-CD38 refractory 116 (89%) 267 (97%) 0.002
Triple refractory 118 (91%) 253 (92%) 0.8
Penta refractory 64 (49%) 141 (51%) 0.7
Prior BCMA-directed therapy 64 (49%) 152 (55%) 0.3
Most recent BCMA

0.012
    CAR-T 54 (84%) 57 (70%)
    TCE 6 (9.4%) 5 (6.2%)
    ADC 4 (6.3%) 19 (23%)
Trial eligible 28 (22%) 63 (23%) 0.8
Received full dose 121 (93%) 268 (97%) 0.093
Baseline ferritin (ng/mL) 293 (132, 832) 393 (131, 1,254) 0.3
Baseline Hgb 10.05 (8.70, 11.70) 9.60 (8.30, 11.30) 0.020
Baseline LDH 223 (177, 277) 208 (166, 287) 0.4
Duration of prior BCMA 0 (0, 22) 0 (0, 126) 0.034
1 Median (Min to Max); n (%); Median (Q1, Q3)
2 Wilcoxon rank sum test; Pearson’s Chi-squared test; Fisher’s exact test

4.3 Unweighted and IPTW baseline characteristics

theme_gtsummary_compact()
## Setting theme "Compact"
data1 = IPTW.data %>%
  transmute(
    "Age (years)" = age_at_bs_ab_initiation,
    "Female sex" = sex_0_male_1_female == 1,
    "Race" = factor(case_when(
      is.na(race_white_0_black_1_aapi_2_american_indian_3_other_4_unk_5) ~ NA,
      race_white_0_black_1_aapi_2_american_indian_3_other_4_unk_5 == 0 ~ "White",
      race_white_0_black_1_aapi_2_american_indian_3_other_4_unk_5 == 1 ~ "Black",
      race_white_0_black_1_aapi_2_american_indian_3_other_4_unk_5 == 2 ~ "AAPI",
      race_white_0_black_1_aapi_2_american_indian_3_other_4_unk_5 == 3 ~ "American Indian",
      race_white_0_black_1_aapi_2_american_indian_3_other_4_unk_5 == 4 ~ "Other",
      race_white_0_black_1_aapi_2_american_indian_3_other_4_unk_5 == 5 ~ NA
    ), levels = c("White","Black","AAPI","American Indian","Other",NA)),
    #ecog_ps_at_bs_ab_initiation_0_5,
    "ECOG" = case_when(
      ecog_ps_at_bs_ab_initiation_0_5 == 0 ~ "0",
      ecog_ps_at_bs_ab_initiation_0_5 == 1 ~ "1",
      ecog_ps_at_bs_ab_initiation_0_5 == 2 ~ "2",
      ecog_ps_at_bs_ab_initiation_0_5 >= 3 ~ "3+"
    ),
    "LDH (U/L)" = as.double(baseline_ldh_prior_to_bs_ab_initiation),
    "Hgb (g/dL)" = baseline_hgb,
    "Baseline ferritin (ng/mL)" = as.double(baseline_ferritin_prior_to_bs_ab_initiation),
    "Heavy chain" = factor(case_when(
      myeloma_involved_heavy_chain_ig_g_ig_a_ig_m_ig_d_ig_e_or_blank == "IgG" ~ "IgG",
      myeloma_involved_heavy_chain_ig_g_ig_a_ig_m_ig_d_ig_e_or_blank == "IgM" ~ "IgM",
      myeloma_involved_heavy_chain_ig_g_ig_a_ig_m_ig_d_ig_e_or_blank == "IgA" ~ "IgA",
      myeloma_involved_heavy_chain_ig_g_ig_a_ig_m_ig_d_ig_e_or_blank == "IgD" ~ "IgD",
      .default = "None"
      ), levels = c("IgG", "IgA","IgM","IgD","None")),
    "True EMD" = just_prior_to_bs_ab_initiation_is_there_true_emd_if_none_move_to_column_gk == 2,
    "≥1 HRCA" = ifelse(deletion_17p_prior_to_bs_ab_0_no_1_yes==1 | t_4_14_prior_to_bs_ab_0_no_1_yes == 1 | t_14_16_prior_to_bs_ab_0_no_1_yes == 1 | amp_1_q_only_0_no_1_yes == 1 | gain_amp1q21_prior_to_bs_ab_0_no_1_yes == 1,1,0),
    "Prior LOTs" = number_of_prior_lines_of_therapy_before_bs_ab,
    "Prior ASCT" = prior_auto_sct_no_0_yes_1,
    "Prior CAR-T" = prior_cart_0_no_1_yes,
    "Triple refractory" = triple_class_refractory_i_mi_d_pi_anti_cd38 == 1,
    "Penta refractory" = penta_refractory_2_p_is_2_i_mi_ds_anti_cd38 == 1,
    "Prior BCMA-directed therapy" = !is.na(most_recent_bcma_agent_name),
    "Trial eligible" = eligible_for_corresponding_bs_ab_trial_on_first_dose_use_trial_elig_sheet_to_help_you == 1,
    "Received full dose" = ifelse(!is.na(date_of_first_bs_ab_full_dose),1,0),
    bs_ab_name,
    weights
  )

table_unweighted <- data1 %>% select(-weights) %>%
  tbl_summary(
    by = bs_ab_name,
    missing = "no",
    #sort = list(all_categorical() ~ "frequency"),
    statistic = list(
      all_continuous() ~ "{median} ({p25}, {p75})",
      all_categorical() ~ "{n} ({p}%)",
      "Age (years)" ~ "{median} ({min} to {max})"
    )
  ) %>%
  add_p() %>%
  bold_labels()

table_IPTW <- survey::svydesign(
    id = ~0,
    weights = ~weights,
    data = data1,
    fpc = NULL
  ) %>%
  tbl_svysummary(
    by = bs_ab_name,
    missing = "no",
    statistic = list(
      all_continuous() ~ "{median} ({p25}, {p75})",
      all_categorical() ~ "{n} ({p}%)"),
    include = c("Age (years)", "Female sex","Race","ECOG","LDH (U/L)","Hgb (g/dL)","Baseline ferritin (ng/mL)", "Heavy chain","True EMD","≥1 HRCA","Prior LOTs","Prior ASCT","Prior CAR-T","Triple refractory","Penta refractory","Trial eligible","Prior BCMA-directed therapy", "Received full dose")
  ) %>%
  add_p() %>%
    bold_labels()
## The following warnings were returned during `add_p()`:
## ! For variable `Heavy chain` (`bs_ab_name`) and "statistic" and "p.value"
##   statistics: Chi-squared approximation may be incorrect
## ! For variable `Race` (`bs_ab_name`) and "statistic" and "p.value" statistics:
##   Chi-squared approximation may be incorrect
tbl_merge(list(table_IPTW, table_unweighted), tab_spanner = c("**IPTW**", "**Unweighted**"))
Characteristic
IPTW
Unweighted
elra
N = 123
1
tec
N = 248
1
p-value2 elra
N = 123
3
tec
N = 248
3
p-value4
Age (years) 70 (64, 76) 69 (62, 76) 0.9 71 (39 to 95) 68 (35 to 88) 0.060
Female sex 66 (54%) 114 (46%) 0.2 68 (55%) 115 (46%) 0.11
Race

0.2

0.15
    White 94 (79%) 187 (76%)
95 (81%) 187 (76%)
    Black 17 (14%) 50 (20%)
16 (14%) 51 (21%)
    AAPI 3 (2.5%) 6 (2.3%)
3 (2.5%) 5 (2.0%)
    American Indian 2 (1.8%) 0 (0%)
2 (1.7%) 0 (0%)
    Other 2 (2.0%) 4 (1.8%)
2 (1.7%) 4 (1.6%)
ECOG

>0.9

0.5
    0 23 (18%) 49 (20%)
21 (17%) 50 (20%)
    1 62 (50%) 119 (48%)
59 (48%) 121 (49%)
    2 28 (23%) 59 (24%)
29 (24%) 60 (24%)
    3+ 10 (8.2%) 20 (8.2%)
14 (11%) 17 (6.9%)
LDH (U/L) 223 (178, 280) 208 (167, 290) 0.4 222 (175, 277) 208 (167, 291) 0.5
Hgb (g/dL) 9.90 (8.50, 11.40) 9.80 (8.40, 11.50) 0.9 10.00 (8.70, 11.70) 9.60 (8.25, 11.20) 0.027
Baseline ferritin (ng/mL) 346 (140, 1,250) 379 (125, 1,219) >0.9 293 (137, 832) 393 (128, 1,254) 0.3
Heavy chain

0.001

<0.001
    IgG 72 (58%) 109 (44%)
72 (59%) 107 (43%)
    IgA 26 (21%) 43 (17%)
26 (21%) 42 (17%)
    IgM 1 (0.7%) 1 (0.4%)
1 (0.8%) 1 (0.4%)
    IgD 2 (1.5%) 0 (0%)
2 (1.6%) 0 (0%)
    None 23 (19%) 95 (38%)
22 (18%) 98 (40%)
True EMD 20 (16%) 39 (16%) >0.9 26 (21%) 34 (14%) 0.067
≥1 HRCA 81 (70%) 151 (65%) 0.3 81 (70%) 152 (65%) 0.4
Prior LOTs 6.00 (5.00, 8.00) 6.00 (4.00, 8.00) 0.4 6.00 (4.00, 8.00) 6.00 (4.00, 8.00) >0.9
Prior ASCT 66 (54%) 141 (57%) 0.6 61 (50%) 144 (58%) 0.12
Prior CAR-T 53 (43%) 90 (36%) 0.2 50 (41%) 94 (38%) 0.6
Triple refractory 112 (91%) 224 (91%) 0.9 112 (91%) 226 (91%) >0.9
Penta refractory 61 (49%) 123 (49%) >0.9 61 (50%) 126 (51%) 0.8
Trial eligible 24 (19%) 55 (22%) 0.5 26 (21%) 53 (21%) >0.9
Prior BCMA-directed therapy 67 (54%) 133 (54%) >0.9 60 (49%) 138 (56%) 0.2
Received full dose 114 (93%) 241 (97%) 0.059 114 (93%) 241 (97%) 0.045
1 Median (Q1, Q3); n (%)
2 Design-based KruskalWallis test; Pearson’s X^2: Rao & Scott adjustment
3 Median (Min to Max); n (%); Median (Q1, Q3)
4 Wilcoxon rank sum test; Pearson’s Chi-squared test; Fisher’s exact test

5 STACKED BAR GRAPH

5.1 FIGURE 1A: Toxicity, IPTW

data1 = IPTW.data %>%
  transmute(
    BsAb = as.factor(bs_ab_name),
    CRS.any = ifelse( max_crs_grade_0_5 != 0, 1, 0),
    ICANS.any = ifelse( max_icans_grade != 0,1,0),
    
    CRS.two = ifelse( max_crs_grade_0_5 == 2 | max_crs_grade_0_5 == 3 | max_crs_grade_0_5 == 4, 1, 0),
    ICANS.two = ifelse( max_icans_grade ==2 | max_icans_grade == 3 | max_icans_grade == 4,1,0),
    
    CRS.three = ifelse( max_crs_grade_0_5 == 3 | max_crs_grade_0_5 == 4, 1, 0),
    ICANS.three = ifelse( max_icans_grade == 3 | max_icans_grade == 4,1,0),
    weights
  )


weighted_table <- survey::svydesign(
    id = ~0,
    weights = ~weights,
    data = data1,
    fpc = NULL
  )


combined_data <- rbind(
  IPTW.data %>%
    count(bs_ab_name, max_crs_grade_0_5, wt=weights) %>%
    group_by(bs_ab_name) %>%
    mutate(
      bs_ab_name = as.factor(recode(bs_ab_name, "tec"="Tec","elra"="Elra")),
      total = sum(n[max_crs_grade_0_5 != 0]),
      pct = prop.table(n) * 100,
      total_pct = sum(pct[max_crs_grade_0_5 != 0]),
      Grade = factor(max_crs_grade_0_5, levels = c(0, 4, 3, 2, 1)),
      Category = "CRS"
    ),
  IPTW.data %>%
    count(bs_ab_name, max_icans_grade, wt=weights) %>%
    group_by(bs_ab_name) %>%
    mutate(
      bs_ab_name = as.factor(recode(bs_ab_name, "tec"="Tec","elra"="Elra")),
      total = sum(n[max_icans_grade != 0]),
      pct = prop.table(n) * 100,
      total_pct = sum(pct[max_icans_grade != 0]),
      Grade = factor(max_icans_grade, levels = c(0, 4, 3, 2, 1)),
      Category = "ICANS"
    )
) %>%
  mutate(
    Category = factor(Category, levels = c("CRS","ICANS"))
  )


# Any grade toxicity
CRS.pval <- as.numeric(svychisq(~CRS.any+BsAb, weighted_table)$p.value)
ICANS.pval <- as.numeric(svychisq(~ICANS.any+BsAb, weighted_table)$p.value)

f_labels <- data.frame(Category = c("CRS","ICANS"), label = c(CRS.pval, ICANS.pval) ) %>%
  mutate(
    Category = factor(Category, levels = c("CRS","ICANS"))
  )

# G2+ toxicity
CRS.two.pval <- as.numeric(svychisq(~CRS.two+BsAb, weighted_table)$p.value)
ICANS.two.pval <- as.numeric(svychisq(~ICANS.two+BsAb, weighted_table)$p.value)

f_labels_two <- data.frame(Category = c("CRS","ICANS"), label = c(CRS.two.pval, ICANS.two.pval) ) %>%
  mutate(
    Category = factor(Category, levels = c("CRS","ICANS"))
  )

# G3+ toxicity
CRS.three.pval <- as.numeric(svychisq(~CRS.three+BsAb, weighted_table)$p.value)
ICANS.three.pval <- as.numeric(svychisq(~ICANS.three+BsAb, weighted_table)$p.value)

f_labels_three <- data.frame(Category = c("CRS","ICANS"), label = c(CRS.three.pval, ICANS.three.pval) ) %>%
  mutate(
    Category = factor(Category, levels = c("CRS","ICANS"))
  )




toxplot_IPTW <- ggplot(combined_data) +
  aes(x = bs_ab_name, y = pct, fill = Grade) +
  geom_bar(stat = "identity", position = "stack", data = . %>% filter(Grade != 0)) +
    geom_text(aes(label = paste0(signif(pct,2), "% (", round(n,0) ,")")), 
            data = . %>% filter(Grade != 0 & Grade != 4),
            position = position_stack(vjust = 0.5)) + 
  geom_text( aes(label = paste0(signif(total_pct,2), "% (", signif(total,2),")"), x = bs_ab_name, y = total_pct + 6, vjust = 0), 
             data = . %>% filter(Grade == 1),
             color = "#104a8e",
             fontface = "bold") +
  facet_wrap(~Category, scales = "free_y") +
  ylab("Proportion of patients (% [n])") +
  xlab(NULL) +
  theme_classic()+
  scale_fill_brewer(palette = "Pastel1") +
  #guides(fill = guide_legend(reverse = TRUE)) +
  scale_y_continuous(breaks = seq(0, 100, by = 20), limits = c(0, 110)) +
  theme(
    plot.title = element_text(size = 14),
    axis.title = element_text(size = 14),
    axis.text = element_text(size = 14),
    legend.title = element_text(size = 14),
    legend.text = element_text(size = 14),
    strip.text = element_text(size = 14),
    axis.text.x = element_text(angle = 45, hjust = 1),
    #legend.position = "bottom"
    legend.position = "right"
  ) + 
  coord_cartesian(ylim = c(0, 115)) +
  geom_text( aes(label = paste0("Grade 1+, p=", signif(label,2)) ), x = 1.5, y = 100+14, vjust = 0, inherit.aes = FALSE, data = f_labels ) + 
  geom_text( aes(label = paste0("Grade 2+, p=", signif(label,2)) ), x = 1.5, y = 100+7, vjust = 0, inherit.aes = FALSE, data = f_labels_two ) + 
  geom_text( aes(label = paste0("Grade 3+, p=", signif(label,2)) ), x = 1.5, y = 100, vjust = 0, inherit.aes = FALSE, data = f_labels_three ) + 
  ggtitle("A) IPTW")

toxplot_IPTW

5.2 FIGURE 1B: Toxicity, unweighted

combined_data <- rbind(
  combined_db %>%
    count(bs_ab_name, max_crs_grade_0_5) %>%
    group_by(bs_ab_name) %>%
    mutate(
      total = sum(n[max_crs_grade_0_5 != 0]),
      pct = prop.table(n) * 100,
      total_pct = sum(pct[max_crs_grade_0_5 != 0]),
      Grade = factor(max_crs_grade_0_5, levels = c(0, 4, 3, 2, 1)),
      bs_ab_name = as.factor(recode(bs_ab_name, "tec"="Tec","elra"="Elra")),
      Category = "CRS"
    ),
  combined_db %>%
    count(bs_ab_name, max_icans_grade) %>%
    group_by(bs_ab_name) %>%
    mutate(
      total = sum(n[max_icans_grade != 0]),
      pct = prop.table(n) * 100,
      total_pct = sum(pct[max_icans_grade != 0]),
      Grade = factor(max_icans_grade, levels = c(0, 4, 3, 2, 1)),
      bs_ab_name = as.factor(recode(bs_ab_name, "tec"="Tec","elra"="Elra")),
      Category = "ICANS"
    )
) %>%
  mutate(
    Category = factor(Category, levels = c("CRS","ICANS"))
    
  )




## Any-grade toxicity
CRS.d.test <- data.frame(
  results = ifelse(combined_db$max_crs_grade_0_5 !=0,1,0) ,
  test = combined_db$bs_ab_name
)
CRS.pval <- chisq.test(CRS.d.test$results,CRS.d.test$test, correct = FALSE)$p.value

ICANS.d.test <- data.frame(
  results = ifelse(combined_db$max_icans_grade !=0,1,0) ,
  test = combined_db$bs_ab_name
)
ICANS.pval <- chisq.test(ICANS.d.test$results,ICANS.d.test$test, correct = FALSE)$p.value


f_labels <- data.frame(Category = c("CRS","ICANS"), label = c(CRS.pval, ICANS.pval) ) %>%
  mutate(
    Category = factor(Category, levels = c("CRS","ICANS") )
  )


## G2+ toxicity
CRS.two.d.test <- data.frame(
  results = ifelse(combined_db$max_crs_grade_0_5 == 2 | combined_db$max_crs_grade_0_5 == 3 | combined_db$max_crs_grade_0_5 == 4,1,0),
  test = combined_db$bs_ab_name
)
CRS.two.pval <- chisq.test(CRS.two.d.test$results,CRS.two.d.test$test, correct = FALSE)$p.value

ICANS.two.d.test <- data.frame(
  results = ifelse(combined_db$max_icans_grade == 2 | combined_db$max_icans_grade == 3 | combined_db$max_icans_grade == 4,1,0),
  test = combined_db$bs_ab_name
)
ICANS.two.pval <- chisq.test(ICANS.two.d.test$results,ICANS.two.d.test$test, correct = FALSE)$p.value

f_labels_two <- data.frame(Category = c("CRS","ICANS"), label = c(CRS.two.pval, ICANS.two.pval) ) %>%
  mutate(
    Category = factor(Category, levels = c("CRS","ICANS") )
  )

## G3+ toxicity
CRS.three.d.test <- data.frame(
  results = ifelse( combined_db$max_crs_grade_0_5 == 3 | combined_db$max_crs_grade_0_5 == 4,1,0),
  test = combined_db$bs_ab_name
)
CRS.three.pval <- fisher.test(CRS.three.d.test$results,CRS.three.d.test$test)$p.value

ICANS.three.d.test <- data.frame(
  results = ifelse(combined_db$max_icans_grade == 3 | combined_db$max_icans_grade == 4,1,0),
  test = combined_db$bs_ab_name
)
ICANS.three.pval <- chisq.test(ICANS.three.d.test$results,ICANS.three.d.test$test, correct = FALSE)$p.value

f_labels_three <- data.frame(Category = c("CRS","ICANS"), label = c(CRS.three.pval, ICANS.three.pval) ) %>%
  mutate(
    Category = factor(Category, levels = c("CRS","ICANS") )
  )



toxplot_unweighted <- ggplot(combined_data) +
  aes(x = bs_ab_name, y = pct, fill = Grade) +
  geom_bar(stat = "identity", position = "stack", data = . %>% filter(Grade != 0)) +
    geom_text(aes(label = paste0(signif(pct,2), "% (",n,")")), 
            data = . %>% filter(Grade != 0 & Grade != 4),
            position = position_stack(vjust = 0.5)) + 
  geom_text( aes(label = paste0(signif(total_pct,2), "% (",total,")"), x = bs_ab_name, y = total_pct + 6, vjust = 0), 
             data = . %>% filter(Grade == 1),
             color = "#104a8e",
             fontface = "bold") +
  facet_wrap(~Category, scales = "free_y") +
  ylab("Proportion of patients (% [n])") +
  xlab(NULL) +
  theme_classic()+
  scale_y_continuous(breaks = seq(0, 100, by = 20), limits = c(0, 100)) +
  theme(
    plot.title = element_text(size = 14),
    axis.title = element_text(size = 14),
    axis.text = element_text(size = 14),
    legend.title = element_text(size = 14),
    legend.text = element_text(size = 14),
    strip.text = element_text(size = 14),
    axis.text.x = element_text(angle = 45, hjust = 1),
    legend.position = "none"
  ) + 
  coord_cartesian(ylim = c(0, 114)) +
  scale_fill_brewer(palette = "Pastel1", direction = 1) +
  geom_text( aes(label = paste0("Grade 1+, p=", signif(label,2)) ), x = 1.5, y = 100+14, vjust = 0, inherit.aes = FALSE, data = f_labels ) + 
  geom_text( aes(label = paste0("Grade 2+, p=", signif(label,2)) ), x = 1.5, y = 100+7, vjust = 0, inherit.aes = FALSE, data = f_labels_two ) + 
  geom_text( aes(label = paste0("Grade 3+, p=", signif(label,2)) ), x = 1.5, y = 100, vjust = 0, inherit.aes = FALSE, data = f_labels_three ) + 
  ggtitle("B) Unweighted")

toxplot_unweighted

5.3 ***FIGURE 1: IPTW and unweighted toicity

FIGURE1 <- toxplot_IPTW / toxplot_unweighted

ggsave("svg/FIGURE1.svg", FIGURE1, device = "svg",
       width = 7,
       height = 11,
       units = c("in"))

ggsave("tiff/FIGURE1.tiff", FIGURE1, device = "tiff",
       width = 7,
       height = 11,
       dpi = 600,
       units = c("in"),
       compression = "lzw")

ggsave("jpeg/FIGURE1.jpeg", FIGURE1, device = "jpeg",
       width = 7,
       height = 11,
       dpi = 600,
       units = "in")

knitr::include_graphics("jpeg/FIGURE1.jpeg")

5.4 FIGURE 2A: Response, IPTW

data1 = IPTW.data %>%
  filter(!is.na(best_response)) %>%
  transmute(
    BsAb = as.factor(bs_ab_name),
    best_response = factor(best_response, levels =  c("sCR", "CR","VGPR","PR","MR","SD","PD")),
    ORR = ifelse(best_response == "sCR" | best_response == "CR" | best_response == "VGPR" | best_response == "PR",1,0),
    VGPR = ifelse(best_response == "sCR" | best_response == "CR" | best_response == "VGPR",1,0),
    CR = ifelse(best_response == "sCR" | best_response == "CR",1,0),
    weights
  )

weighted_table <- survey::svydesign(
    id = ~0,
    weights = ~weights,
    data = data1,
    fpc = NULL
  )



data1 <-   IPTW.data %>%
  filter(!is.na(best_response)) %>%
  count(bs_ab_name, best_response, wt=weights) %>%
  group_by(bs_ab_name) %>%
  mutate(
    total = sum(n[best_response == "sCR" | best_response == "CR" | best_response == "VGPR" | best_response == "PR"]),
    pct = prop.table(n) * 100,
    total_pct = sum(pct[best_response == "sCR" | best_response == "CR" | best_response == "VGPR" | best_response == "PR"]),
    best_response = factor(best_response, levels =  c("sCR", "CR","VGPR","PR","MR", "SD","PD")),
    bs_ab_name = as.factor(bs_ab_name),
    Category = "Response rate"
  )



Response.pval <- as.numeric(svychisq(~ORR+BsAb, weighted_table)$p.value)
CR.pval <- as.numeric(svychisq(~CR+BsAb, weighted_table)$p.value)
VGPR.pval <- as.numeric(svychisq(~VGPR+BsAb, weighted_table)$p.value)

f_labels <- data.frame(Category = c("Response rate"), label = c(Response.pval) ) %>%
  mutate(
    Category = factor(Category, levels = c("Response rate"))
  )

f_labels_CR <- data.frame(Category = c("Response rate"), label = c(CR.pval) ) %>%
  mutate(
    Category = factor(Category, levels = c("Response rate"))
  )

f_labels_VGPR <- data.frame(Category = c("Response rate"), label = c(VGPR.pval) ) %>%
  mutate(
    Category = factor(Category, levels = c("Response rate"))
  )

responseplot_IPTW <- ggplot(data1) +
  aes(x = bs_ab_name, y = pct, fill = best_response) +
  geom_bar(stat = "identity", position = "stack", data = . %>% filter(best_response == "sCR" | best_response == "CR" | best_response == "VGPR" | best_response == "PR")) +
    geom_text(aes(label = paste0(signif(pct,2), "% (", round(n,0) ,")")), 
            data = . %>% filter(best_response == "sCR" | best_response == "CR" | best_response == "VGPR"| best_response == "PR"),
            position = position_stack(vjust = 0.5)) + 
  geom_text( aes(label = paste0(signif(total_pct,2), "% (", signif(total,2),")"), x = bs_ab_name, y = total_pct + 7, vjust = 0), 
             data = . %>% filter(best_response == "PR"),
             color = "#104a8e",
             fontface = "bold") +
  facet_wrap(~Category, scales = "free_y") +
  ylab("Proportion of patients (% [n])") +
  xlab(NULL) +
  theme_classic()+
  scale_fill_brewer(
    palette = "Pastel1",
    direction = -1,
    breaks = rev(levels(data1$best_response))
  ) +
  guides(fill = guide_legend(reverse = TRUE)) +
  scale_y_continuous(breaks = seq(0, 100, by = 20), limits = c(0, 110)) +
  theme(
    plot.title = element_text(size = 14),
    axis.title = element_text(size = 14),
    axis.text = element_text(size = 14),
    legend.title = element_text(size = 14),
    legend.text = element_text(size = 14),
    strip.text = element_text(size = 14),
    axis.text.x = element_text(angle = 45, hjust = 1),
    #legend.position = "bottom"
    legend.position = "right"
  ) + 
  coord_cartesian(ylim = c(0, 114)) +
  geom_text( aes(label = paste0("ORR, p=", signif(Response.pval,2)) ), x = 1.5, y = 100+14, vjust = 0, inherit.aes = FALSE, data = f_labels  ) + 
  geom_text( aes(label = paste0("≥VGPR, p=", signif(VGPR.pval,2)) ), x = 1.5, y = 100+7, vjust = 0, inherit.aes = FALSE, data = f_labels_VGPR  ) + 
  geom_text( aes(label = paste0("≥CR, p=", signif(CR.pval,2)) ), x = 1.5, y = 100, vjust = 0, inherit.aes = FALSE, data = f_labels_CR  ) + 
  labs(fill = "Best response")+ 
  ggtitle("A) IPTW")

responseplot_IPTW

5.5 FIGURE 2B: Response, unweighted

data1 <- combined_db %>%
  filter(!is.na(best_response)) %>%
  transmute(
    best_response = factor(
      best_response, 
      levels =  c("sCR", "CR","VGPR", "PR","MR","SD","PD")),
    ORR = ifelse(best_response == "sCR" | best_response == "CR" | best_response == "VGPR" |best_response == "PR",1,0),
    bs_ab_name = as.factor(recode(bs_ab_name, "tec"="Tec","elra"="Elra")),
    Category = "Response rate"
  )

combined_data <- combined_db %>%
  filter(!is.na(best_response)) %>%
  count(bs_ab_name, best_response) %>%
  group_by(bs_ab_name) %>%
  transmute(
    bs_ab_name = as.factor(recode(bs_ab_name, "tec"="Tec","elra"="Elra")),
    best_response = factor(
      best_response, 
      levels =  c("sCR", "CR","VGPR", "PR","MR","SD","PD")),
    total = sum(n[best_response == "sCR" | best_response == "CR" | best_response == "VGPR" | best_response == "PR"]),
    pct = prop.table(n) * 100,
    n,
    total_pct = sum(pct[best_response == "sCR" | best_response == "CR" |best_response == "VGPR" | best_response == "PR"]),
    Category = "Response rate"
  )

## ORR
Response.d.test <- data.frame(
  results = ifelse(data1$ORR,1,0),
  test = data1$bs_ab_name
)
Response.pval <- chisq.test(Response.d.test$results,Response.d.test$test, correct = FALSE)$p.value

## VGPR
VGPR.d.test <- data.frame(
  results = ifelse(data1$best_response == "CR" | data1$best_response == "VGPR",1,0),
  test = data1$bs_ab_name
)
VGPR.pval <- chisq.test(VGPR.d.test$results,VGPR.d.test$test, correct = FALSE)$p.value

## CR
CR.d.test <- data.frame(
  results = ifelse(data1$best_response == "CR",1,0),
  test = data1$bs_ab_name
)
CR.pval <- chisq.test(CR.d.test$results,CR.d.test$test, correct = FALSE)$p.value


f_labels <- data.frame(Category = c("Response rate"), label = c(Response.pval) ) %>%
  mutate(
    Category = factor(Category, levels = c("Response rate"))
  )

f_labels_CR <- data.frame(Category = c("Response rate"), label = c(CR.pval) ) %>%
  mutate(
    Category = factor(Category, levels = c("Response rate"))
  )

f_labels_VGPR <- data.frame(Category = c("Response rate"), label = c(VGPR.pval) ) %>%
  mutate(
    Category = factor(Category, levels = c("Response rate"))
  )


responseplot_unweighted <- ggplot(combined_data) +
  aes(x = bs_ab_name, y = pct, fill = best_response) +
  geom_bar(stat = "identity", position = "stack", data = . %>% filter(best_response == "sCR" | best_response == "CR" |best_response == "VGPR" | best_response == "PR")) +
  geom_text(aes(label = paste0(signif(pct,2), "% (",n,")")), 
             data = . %>% filter(best_response == "sCR" | best_response == "CR" | best_response == "VGPR" | best_response == "PR"),
             position = position_stack(vjust = 0.5)) + 
  geom_text( aes(label = paste0(signif(total_pct,2), "% (",total,")"), x = bs_ab_name, y = total_pct + 7, vjust = 0), 
             data = . %>% filter(best_response == "PR"),
             color = "#104a8e",
             fontface = "bold") +
  facet_wrap(~Category, scales = "free_y") +
  ylab("Proportion of patients (% [n])") +
  xlab(NULL) +
  theme_classic()+
  scale_y_continuous(breaks = seq(0, 100, by = 20), limits = c(0, 110)) +
  theme(
    plot.title = element_text(size = 14),
    axis.title = element_text(size = 14),
    axis.text = element_text(size = 14),
    legend.title = element_text(size = 14),
    legend.text = element_text(size = 14),
    strip.text = element_text(size = 14),
    axis.text.x = element_text(angle = 45, hjust = 1),
    #legend.position = "bottom"
    legend.position = "none"
  ) + 
  coord_cartesian(ylim = c(0, 114)) +
  scale_fill_brewer(palette = "Pastel1", direction = -1) +
  geom_text( aes(label = paste0("ORR, p=", signif(Response.pval,2)) ), x = 1.5, y = 100+14, vjust = 0, inherit.aes = FALSE, data = f_labels) + 
  geom_text( aes(label = paste0("≥VGPR, p=", signif(VGPR.pval,2)) ), x = 1.5, y = 100+7, vjust = 0, inherit.aes = FALSE, data = f_labels_VGPR  ) + 
  geom_text( aes(label = paste0("≥CR, p=", signif(CR.pval,2)) ), x = 1.5, y = 100, vjust = 0, inherit.aes = FALSE, data = f_labels_CR  ) + 
  labs(fill = "Best response")+ 
  ggtitle("B) Unweighted")

responseplot_unweighted

5.6 ***FIGURE 2: IPTW and unweighted response

FIGURE2 <- responseplot_IPTW / responseplot_unweighted

ggsave("svg/FIGURE2.svg", FIGURE2, device = "svg",
       width = 7,
       height = 11,
       units = c("in"))

ggsave("tiff/FIGURE2.tiff", FIGURE2, device = "tiff",
       width = 7,
       height = 11,
       dpi = 600,
       units = c("in"),
       compression = "lzw")

ggsave("jpeg/FIGURE2.jpeg", FIGURE2, device = "jpeg",
       width = 7,
       height = 11,
       dpi = 600,
       units = "in")

knitr::include_graphics("jpeg/FIGURE2.jpeg")

6 SURVIVAL ANALYSIS

6.1 Length of follow-up (Reverse KM, months)

data1 <- IPTW.data %>%
  transmute(
    bs_ab_name,
    Death = ifelse(death_yes_no,1,0),
    Days.to.death.or.DLC = as.numeric(date_of_last_contact - date_of_bs_ab_step_up_d1),
    Reverse_death = ifelse(death_yes_no == 1, 0,1),
    weights
  )

reverse_km_OS <- survfit(Surv(Days.to.death.or.DLC/30, Reverse_death) ~ 1, weights = weights, data1)
reverse_km_OS
## Call: survfit(formula = Surv(Days.to.death.or.DLC/30, Reverse_death) ~ 
##     1, data = data1, weights = weights)
## 
##    15 observations deleted due to missingness 
##      records   n events median 0.95LCL 0.95UCL
## [1,]     356 356    215   13.4    12.5    14.9
reverse_km_OS_elra <- survfit(Surv(Days.to.death.or.DLC/30, Reverse_death) ~ 1, weights = weights, data1 %>% filter(bs_ab_name == "elra"))
reverse_km_OS_elra
## Call: survfit(formula = Surv(Days.to.death.or.DLC/30, Reverse_death) ~ 
##     1, data = data1 %>% filter(bs_ab_name == "elra"), weights = weights)
## 
##    6 observations deleted due to missingness 
##      records   n events median 0.95LCL 0.95UCL
## [1,]     117 117   74.5   7.47    6.03       9
reverse_km_OS_tec <- survfit(Surv(Days.to.death.or.DLC/30, Reverse_death) ~ 1, weights = weights, data1 %>% filter(bs_ab_name == "tec"))
reverse_km_OS_tec
## Call: survfit(formula = Surv(Days.to.death.or.DLC/30, Reverse_death) ~ 
##     1, data = data1 %>% filter(bs_ab_name == "tec"), weights = weights)
## 
##    9 observations deleted due to missingness 
##      records   n events median 0.95LCL 0.95UCL
## [1,]     239 239    140   17.1    14.9    20.1

6.2 Stratified by BsAb

data1 <- IPTW.data %>%
  transmute(
    site,
    patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate,
    bs_ab_name,
    Relapse.or.death = ifelse(progression_yes_no | death_yes_no,1,0),
    Death = ifelse(death_yes_no,1,0),
    date_of_bs_ab_step_up_d1,
    date_of_pd,
    date_of_last_contact,
    date_progression_or_death = coalesce(date_of_pd, date_of_last_contact),
    Days.to.death.or.DLC = as.numeric(date_of_last_contact - date_of_bs_ab_step_up_d1),
    Days.to.death.or.relapse.or.DLC = as.numeric(date_progression_or_death  - date_of_bs_ab_step_up_d1),
    trial_eligible = ifelse(eligible_for_corresponding_bs_ab_trial_on_first_dose_use_trial_elig_sheet_to_help_you ==1,1,0),
    weights
  )

data_dor <- IPTW.data %>%
  filter(best_response %in% c("sCR","CR","VGPR","PR")) %>%
  transmute(
    weights,
    bs_ab_name,
    day_30_response,
    day_90_response,
    best_response,
    Relapse.or.death = ifelse(progression_yes_no | death_yes_no,1,0),
    Death = ifelse(death_yes_no,1,0),
    date_of_pd,
    date_of_last_contact,
    date_of_bs_ab_step_up_d1,
    date_progression_or_death = coalesce(date_of_pd, date_of_last_contact),
    day_30_response,
    day_90_response,
    best_response,
    time_to_best_response = case_when(
      best_response == day_30_response ~ 30,
      best_response == day_90_response ~ 90,
      .default = NA
    ),
    trial_eligible = ifelse(eligible_for_corresponding_bs_ab_trial_on_first_dose_use_trial_elig_sheet_to_help_you ==1,1,0)
  ) %>%
  filter(!is.na(time_to_best_response)) %>%
  mutate(
    Days.to.death.or.DLC = as.numeric(date_of_last_contact - date_of_bs_ab_step_up_d1) - time_to_best_response,
    Days.to.death.or.relapse.or.DLC = as.numeric(date_progression_or_death  - date_of_bs_ab_step_up_d1)- time_to_best_response,
  ) %>%
  filter(Days.to.death.or.DLC >0 & Days.to.death.or.relapse.or.DLC >0)

km_OS <- survfit(Surv(Days.to.death.or.DLC/30, Death) ~ bs_ab_name, weights = weights, data = data1)
km_OS
## Call: survfit(formula = Surv(Days.to.death.or.DLC/30, Death) ~ bs_ab_name, 
##     data = data1, weights = weights)
## 
##    15 observations deleted due to missingness 
##                 records   n events median 0.95LCL 0.95UCL
## bs_ab_name=elra     117 117   42.8     NA     8.6      NA
## bs_ab_name=tec      239 239   98.6   21.5    16.4      NA
med_time <- surv_median(km_OS)

cox_model <- coxph(Surv(Days.to.death.or.DLC/30, Death) ~ ifelse(bs_ab_name == "elra",1,0), weights = weights, data = data1)
cox_summary <- summary(cox_model)

HR <- round(cox_summary$coefficients[1, "exp(coef)"], 2)
CI_lower <- round(cox_summary$conf.int[1, "lower .95"], 2)
CI_upper <- round(cox_summary$conf.int[1, "upper .95"], 2)
pval <- signif(cox_summary$coefficients[1, "Pr(>|z|)"], 2)

## 1-year timepoint
timepoint1 <- summary(km_OS,times=c(3))
timepoint2 <- summary(km_OS,times=c(6))
timepoint3 <- summary(km_OS,times=c(9))
timepoint4 <- summary(km_OS,times=c(12))
timepoint5 <- summary(km_OS,times=c(15))

OS <- ggsurvplot(km_OS,
           pval=FALSE,
           conf.int = FALSE,
           risk.table = TRUE, # Add risk table
           fontsize = 4,
           risk.table.col = "strata", # Change risk table color by groups
           tables.height = 0.3,
           tables.theme = theme_cleantable() +
             theme(
               axis.text.y = element_text(size = 9),
               plot.title = element_text(size = 10)
               ),
           xlab="Time (months)", ylab = "Overall survival",
           xlim = c(0,18), break.x.by = c(3),
           ylim = c(0,1), break.y.by = c(0.25),
           legend = "none",
           surv.scale="percent",
           font.main = c(10, "plain", "black"),
           font.x = c(10, "plain", "black"),
           font.y = c(10, "plain", "black"),
           font.caption = c(10, "plain", "black"),
           font.tickslab = c(10, "plain", "black"),
           legend.labs = c("Elra", "Tec"),
           title = "A)",
           size = 0.5
) 

OS$plot <- OS$plot + 
  geom_vline(xintercept = c(3, 6, 9, 12, 15), linetype = "dotted", color = "gray50") +
  annotate("text", x = 3, y = timepoint1$surv[1] - 0.1, label = paste0(signif(timepoint1$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 3, y = timepoint1$surv[2] + 0.03, label = paste0(signif(timepoint1$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 6, y = timepoint2$surv[1] - 0.08, label = paste0(signif(timepoint2$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 6, y = timepoint2$surv[2] + 0.03, label = paste0(signif(timepoint2$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 9, y = timepoint3$surv[1] - 0.08, label = paste0(signif(timepoint3$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 9, y = timepoint3$surv[2] + 0.03, label = paste0(signif(timepoint3$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 12, y = timepoint4$surv[1] - 0.08, label = paste0(signif(timepoint4$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 12, y = timepoint4$surv[2] + 0.03, label = paste0(signif(timepoint4$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 15, y = timepoint5$surv[1] - 0.08, label = paste0(signif(timepoint5$surv[1]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 15, y = timepoint5$surv[2] + 0.03, label = paste0(signif(timepoint5$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = -0.7, y = 0,
           label = paste0("Elra vs tec\n",
                          " Median: ", signif(med_time$median[1], 3), " vs ", signif(med_time$median[2], 3), " months\n",
                          " HR: ", sprintf("%.2f", HR), " (95% CI ", sprintf("%.2f", CI_lower), "-", sprintf("%.2f", CI_upper),", p=",pval, ")"),
           color = "black", hjust = 0, vjust = 0, size = 3)

OS

km_PFS <- survfit(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ bs_ab_name, weights = weights, data = data1)
km_PFS
## Call: survfit(formula = Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ 
##     bs_ab_name, data = data1, weights = weights)
## 
##    12 observations deleted due to missingness 
##                 records   n events median 0.95LCL 0.95UCL
## bs_ab_name=elra     118 118   68.1   3.77    2.57    7.17
## bs_ab_name=tec      241 241  151.1   8.50    7.10   11.30
med_time <- surv_median(km_PFS)

cox_model <- coxph(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ ifelse(bs_ab_name == "elra",1,0), weights = weights, data = data1)
cox_summary <- summary(cox_model)

HR <- round(cox_summary$coefficients[1, "exp(coef)"], 2)
CI_lower <- round(cox_summary$conf.int[1, "lower .95"], 2)
CI_upper <- round(cox_summary$conf.int[1, "upper .95"], 2)
pval <- signif(cox_summary$coefficients[1, "Pr(>|z|)"], 2)

## Timepoints
timepoint1 <- summary(km_PFS,times=c(3))
timepoint2 <- summary(km_PFS,times=c(6))
timepoint3 <- summary(km_PFS,times=c(9))
timepoint4 <- summary(km_PFS,times=c(12))
timepoint5 <- summary(km_PFS,times=c(15))


PFS <- ggsurvplot(km_PFS,
           pval=FALSE,
           conf.int = FALSE,
           risk.table = TRUE, # Add risk table
           fontsize = 4,
           risk.table.col = "strata", # Change risk table color by groups
           tables.height = 0.3,
           tables.theme = theme_cleantable() +
             theme(
               axis.text.y = element_text(size = 9),
               plot.title = element_text(size = 10)
               ),
           xlab="Time (months)", ylab = "Progression-free survival",
           xlim = c(0,18), break.x.by = c(3),
           ylim = c(0,1), break.y.by = c(0.25),
           legend = "none",
           surv.scale="percent",
           font.main = c(10, "plain", "black"),
           font.x = c(10, "plain", "black"),
           font.y = c(10, "plain", "black"),
           font.caption = c(10, "plain", "black"),
           font.tickslab = c(10, "plain", "black"),
           legend.labs = c("Elra", "Tec"),
           title = "B)",
           size = 0.5
)

PFS$plot <- PFS$plot + 
  geom_vline(xintercept = c(3, 6, 9, 12, 15), linetype = "dotted", color = "gray50") +
  annotate("text", x = 3, y = timepoint1$surv[1] - 0.08, label = paste0(signif(timepoint1$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 3, y = timepoint1$surv[2] + 0.03, label = paste0(signif(timepoint1$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 6, y = timepoint2$surv[1] - 0.08, label = paste0(signif(timepoint2$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 6, y = timepoint2$surv[2] + 0.03, label = paste0(signif(timepoint2$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 9, y = timepoint3$surv[1] - 0.08, label = paste0(signif(timepoint3$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 9, y = timepoint3$surv[2] + 0.03, label = paste0(signif(timepoint3$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 12, y = timepoint4$surv[1] - 0.08, label = paste0(signif(timepoint4$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 12, y = timepoint4$surv[2] + 0.03, label = paste0(signif(timepoint4$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 15, y = timepoint5$surv[1] - 0.08, label = paste0(signif(timepoint5$surv[1]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 15, y = timepoint5$surv[2] + 0.03, label = paste0(signif(timepoint5$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = -0.7, y = 0,
           label = paste0("Elra vs tec\n",
                          " Median: ", signif(med_time$median[1], 3), " vs ", signif(med_time$median[2], 3), " months\n",
                          " HR: ", sprintf("%.2f", HR), " (95% CI ", sprintf("%.2f", CI_lower), "-", sprintf("%.2f", CI_upper),", p=",pval, ")"),
           color = "black", hjust = 0, vjust = 0, size = 3)

PFS

km_DOR <- survfit(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ bs_ab_name, weights = weights, data = data_dor)
km_DOR
## Call: survfit(formula = Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ 
##     bs_ab_name, data = data_dor, weights = weights)
## 
##    3 observations deleted due to missingness 
##                 records     n events median 0.95LCL 0.95UCL
## bs_ab_name=elra      62  60.2   20.6   11.7    6.13      NA
## bs_ab_name=tec      118 120.7   56.1   13.7    9.27    22.3
med_time <- surv_median(km_DOR)

cox_model <- coxph(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ ifelse(bs_ab_name == "elra",1,0), weights = weights, data = data_dor)
cox_summary <- summary(cox_model)

HR <- round(cox_summary$coefficients[1, "exp(coef)"], 2)
CI_lower <- round(cox_summary$conf.int[1, "lower .95"], 2)
CI_upper <- round(cox_summary$conf.int[1, "upper .95"], 2)
pval <- signif(cox_summary$coefficients[1, "Pr(>|z|)"], 2)

## Specific timepoints
timepoint1 <- summary(km_DOR,times=c(3))
timepoint2 <- summary(km_DOR,times=c(6))
timepoint3 <- summary(km_DOR,times=c(9))
timepoint4 <- summary(km_DOR,times=c(12))
timepoint5 <- summary(km_DOR,times=c(15))


DOR <- ggsurvplot(km_DOR,
           pval=FALSE,
           conf.int = FALSE,
           risk.table = TRUE, # Add risk table
           fontsize = 4,
           risk.table.col = "strata", # Change risk table color by groups
           tables.height = 0.3,
           #risk.table.fontsize = 4,
           tables.theme = theme_cleantable() +
             theme(
               axis.text.y = element_text(size = 9),
               plot.title = element_text(size = 10)
               ),
           xlab="Time (months)", ylab = "Duration of response",
           xlim = c(0,18), break.x.by = c(3),
           ylim = c(0,1), break.y.by = c(0.25),
           legend = "none",
           surv.scale="percent",
           font.main = c(10, "plain", "black"),
           font.x = c(10, "plain", "black"),
           font.y = c(10, "plain", "black"),
           font.caption = c(10, "plain", "black"),
           font.tickslab = c(10, "plain", "black"),
           legend.labs = c("Elra", "Tec"),
           title = "C)",
           size = 0.5
)

DOR$plot <- DOR$plot + 
  geom_vline(xintercept = c(3, 6, 9, 12, 15), linetype = "dotted", color = "gray50") +
  annotate("text", x = 3, y = timepoint1$surv[1] -0.08, label = paste0(signif(timepoint1$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 3, y = timepoint1$surv[2] + 0.03, label = paste0(signif(timepoint1$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 6, y = timepoint2$surv[1] - 0.08, label = paste0(signif(timepoint2$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 6, y = timepoint2$surv[2] + 0.03, label = paste0(signif(timepoint2$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 9, y = timepoint3$surv[1] - 0.08, label = paste0(signif(timepoint3$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 9, y = timepoint3$surv[2] + 0.03, label = paste0(signif(timepoint3$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 12, y = timepoint4$surv[1] - 0.08, label = paste0(signif(timepoint4$surv[1]*100, 2), "%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 12, y = timepoint4$surv[2] + 0.03, label = paste0(signif(timepoint4$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 15, y = timepoint5$surv[1] - 0.08, label = paste0(signif(timepoint5$surv[1]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = 15, y = timepoint5$surv[2] + 0.03, label = paste0(signif(timepoint5$surv[2]*100, 2),"%"), color = "black", hjust = 0, vjust = 0, size = 3) +
  annotate("text", x = -0.7, y = 0,
           label = paste0("Elra vs tec\n",
                          " Median: ", signif(med_time$median[1], 3), " vs ", signif(med_time$median[2], 3), " months\n",
                          " HR: ", sprintf("%.2f", HR), " (95% CI ", sprintf("%.2f", CI_lower), "-", sprintf("%.2f", CI_upper),", p=",pval, ")"),
           color = "black", hjust = 0, vjust = 0, size = 3)

DOR

6.3 ***FIGURE 3: PFS/OS/DOR, stratified by trial elibility

FIGURE3 <- wrap_plots(
  (OS$plot / OS$table) + plot_layout(heights = c(4, 0.9)),
  (PFS$plot / PFS$table) + plot_layout(heights = c(4, 0.9)),
  (DOR$plot / DOR$table) + plot_layout(heights = c(4, 0.9)),
  ncol = 1  # Force vertical stacking
)

ggsave("svg/FIGURE3.svg", FIGURE3, device = "svg",
       width = 6.5,
       height = 8.5,
       units = c("in"))

ggsave("jpeg/FIGURE3.jpeg", FIGURE3, device = "jpeg",
       width = 6.5,
       height = 8.5,
       dpi = 600,
       units = c("in"))

ggsave("tiff/FIGURE3.tiff", FIGURE3, device = "tiff",
       width = 6.5,
       height = 8.5,
       dpi = 600,
       units = c("in"),
       compression = "lzw")

knitr::include_graphics("jpeg/FIGURE3.jpeg")

6.4 Schoenfeld residuals and PH test

data1 <- IPTW.data %>%
  transmute(
    site,
    patient_identifier_use_same_id_for_cns_amyloid_sheets_if_appropriate,
    elra_ind = ifelse(bs_ab_name == "elra", 1, 0),  # 1 = elra, 0 = tec
    Relapse.or.death = ifelse(progression_yes_no | death_yes_no,1,0),
    Death = ifelse(death_yes_no,1,0),
    date_of_bs_ab_step_up_d1,
    date_of_pd,
    date_of_last_contact,
    date_progression_or_death = coalesce(date_of_pd, date_of_last_contact),
    Days.to.death.or.DLC = as.numeric(date_of_last_contact - date_of_bs_ab_step_up_d1),
    Days.to.death.or.relapse.or.DLC = as.numeric(date_progression_or_death  - date_of_bs_ab_step_up_d1),
    trial_eligible = ifelse(eligible_for_corresponding_bs_ab_trial_on_first_dose_use_trial_elig_sheet_to_help_you ==1,1,0),
    weights
  )

## ---- Cox models ----
cox_OS <- coxph(Surv(Days.to.death.or.DLC/30, Death) ~ elra_ind, weights = weights, data = data1)

cox_PFS <- coxph(Surv(Days.to.death.or.relapse.or.DLC/30, Relapse.or.death) ~ elra_ind, weights = weights, data = data1)

## ---- PH diagnostics ----
zph_OS  <- cox.zph(cox_OS)
zph_PFS <- cox.zph(cox_PFS)

# Global PH p-values (same as elra_ind right now)
p_OS  <- signif(zph_OS$table["GLOBAL",  "p"], 3)
p_PFS <- signif(zph_PFS$table["GLOBAL", "p"], 3)

par(mfrow = c(1, 2))

plot(
  zph_OS,
  main = paste0(
    "Schoenfeld residuals – OS\n",
    "Global PH test: p = ", p_OS
  )
)

plot(
  zph_PFS,
  main = paste0(
    "Schoenfeld residuals – PFS\n",
    "Global PH test: p = ", p_PFS
  )
)

par(mfrow = c(1, 1))

7 REGRESSION

7.1 ORR, ≥CR, ≥VGPR, PFS, OS, DOR, IFS, CRS G1+, ICANS G1+, steroid use (IPTW)

options(gtsummary.add_p_footnote = FALSE, gtsummary.add_estimate_to_reference_rows = FALSE)
theme_gtsummary_compact()
## Setting theme "Compact"
preds0 <- c("Response", "BsAb", "Baseline Hgb", "Baseline LDH", "ECOG 2+", "True EMD", "Prior LOTs", "Prior BCMA")
preds1 <- c("CR/sCR", "BsAb", "Baseline Hgb", "Baseline LDH", "ECOG 2+", "True EMD", "Prior LOTs", "Prior BCMA")
preds2 <- c("BsAb", "Baseline Hgb",  "Baseline LDH", "ECOG 2+","True EMD", "Prior LOTs", "Prior BCMA")
preds3 <- c("ICANS 1+", "BsAb", "Baseline Hgb", "Baseline LDH", "ECOG 2+", "True EMD", "Prior LOTs", "Prior BCMA")
preds4 <- c("CRS 1+", "BsAb", "Baseline Hgb", "Baseline LDH", "ECOG 2+", "True EMD", "Prior LOTs", "Prior BCMA")
preds5 <- c("Steroid use", "BsAb", "Baseline Hgb", "Baseline LDH", "ECOG 2+", "True EMD", "Prior LOTs", "Prior BCMA")

preds6 <- c("VGPR", "BsAb", "Baseline Hgb", "Baseline LDH", "ECOG 2+", "True EMD", "Prior LOTs", "Prior BCMA")

data0 <- IPTW.data %>%
  transmute(
    weights,
    BsAb = factor(recode(bs_ab_name, "tec"="Tec","elra"="Elra"), levels = c("Tec","Elra")),
    HT_Risk,
    best_response,
    "Infection" = ifelse(infection == 1, 1,0),
    "CrCl <40" = ifelse(cr_cl_prior_to_step_up_dose_number_1 <40, 1,0),
    "CR/sCR" = ifelse(best_response %in% c("CR","sCR"),1,0),
    "VGPR" = ifelse(best_response %in% c("CR","sCR", "VGPR"),1,0),
    "Response" = ifelse(best_response %in% c("CR","sCR", "VGPR", "PR"),1,0),
    "True EMD" = ifelse(just_prior_to_bs_ab_initiation_is_there_true_emd_if_none_move_to_column_gk == 2, 1,0),
    "Steroid use" = ifelse(steroid_use_as_tx==1,1,0),
    "Age" = age_at_bs_ab_initiation/10,
    "ECOG" = ecog_ps_at_bs_ab_initiation_0_5,
    "ECOG 2+" = case_when(
      is.na(ecog_ps_at_bs_ab_initiation_0_5) ~ NA,
      ecog_ps_at_bs_ab_initiation_0_5 %in% c(2,3,4) ~ 1,
      TRUE ~ 0),
    "ICANS 1+" = case_when(
      is.na(max_icans_grade) ~ NA,
      max_icans_grade %in% c(1,2,3,4) ~1,
      TRUE ~ 0
      ),
    "CRS 1+" = case_when(
      is.na(max_crs_grade_0_5) ~ NA,
      max_crs_grade_0_5 %in% c(1,2,3,4) ~1,
      TRUE ~ 0
      ),
    "Prior LOTs" = number_of_prior_lines_of_therapy_before_bs_ab,
    "HRCA" = ifelse(deletion_17p_prior_to_bs_ab_0_no_1_yes | t_4_14_prior_to_bs_ab_0_no_1_yes | t_14_16_prior_to_bs_ab_0_no_1_yes | gain_amp1q21_prior_to_bs_ab_0_no_1_yes, 1, 0),
    "Prior BCMA" = prior_bcma_directed_therapy_yes_no,
    prior_bcma_directed_therapy_yes_no,
    Days.since.last.BCMA = case_when(
      is.na(prior_bcma_directed_therapy_yes_no) ~ NA,
      date_of_last_exposure_to_prior_bcma_agent > date_of_bs_ab_step_up_d1 ~ NA,
      prior_bcma_directed_therapy_yes_no == 1 ~ as.numeric(as.Date(date_of_bs_ab_step_up_d1) - as.Date(date_of_last_exposure_to_prior_bcma_agent))
    ),
    Relapse.or.death = ifelse(progression_yes_no | death_yes_no,1,0),
    Death = ifelse(death_yes_no,1,0),
    date_of_pd,
    date_of_last_contact,
    date_progression_or_death = coalesce(date_of_pd, date_of_last_contact),
    Days.to.death.or.DLC = as.numeric(date_of_last_contact - date_of_bs_ab_step_up_d1),
    Days.to.death.or.relapse.or.DLC = as.numeric(date_progression_or_death  - date_of_bs_ab_step_up_d1),
    "ALC BL" = baseline_alc,
    "CrCl BL" = cr_cl_prior_to_step_up_dose_number_1,
    "Trial ineligible" = ifelse(eligible_for_corresponding_bs_ab_trial_on_first_dose_use_trial_elig_sheet_to_help_you == 1,0,1),
    "Baseline platelets" = baseline_pl_ts,
    "Baseline Hgb" = baseline_hgb,
    "Baseline LDH" = as.double(baseline_ldh_prior_to_bs_ab_initiation)/100,
    "Baseline ferritin" = as.double(baseline_ldh_prior_to_bs_ab_initiation)/1000,
    time_to_best_response = case_when(
      best_response == day_30_response ~ 30,
      best_response == day_90_response ~ 90,
      .default = NA
      ),
    date_of_bs_ab_step_up_d1,
    bs_ab_served_only_as_pre_car_t_bridging_0_no_1_yes,
    days_to_infection_or_last_contact = as.numeric(
      pmin(as.Date(date_of_infection_74), as.Date(date_of_last_contact), na.rm = TRUE) - as.Date(date_of_bs_ab_step_up_d1)
    ),
    infection_or_death = ifelse(!is.na(date_of_infection_74) | death_yes_no == 1,1,0),
    IVIg = ifelse(ivig_use_yes_1_no_0 == 1, 1, 0),
    ALPS = ifelse(baseline_hgb >10 & baseline_ldh_prior_to_bs_ab_initiation <300,1,0),
    baseline_pl_ts,
    baseline_anc,
    baseline_hgb,
    baseline_crp_prior_to_bs_ab_initiation,
    baseline_ferritin_prior_to_bs_ab_initiation,
    CAR.HEMATOTOX = baseline_pl_ts * -0.05 + baseline_anc * -0.03 + baseline_hgb * -0.20 + as.numeric(baseline_crp_prior_to_bs_ab_initiation) * 0.15 + as.numeric(baseline_ferritin_prior_to_bs_ab_initiation) * 0.001
    ) %>%
  mutate(
    "BCMA exposure status" = factor(case_when(
      #prior_bcma_directed_therapy_yes_no == 0 ~ "  Never",
      Days.since.last.BCMA < 365 ~ "  <1 yr",
      Days.since.last.BCMA >= 365 ~ "  ≥1 yr"
    ), c("  ≥1 yr","  <1 yr")),
    `CAR-HEMATOTOX risk` = factor(
      case_when(
        CAR.HEMATOTOX < 0 ~ "Low",
        CAR.HEMATOTOX <= 5 ~ "Intermediate",
        CAR.HEMATOTOX > 5 ~ "High"
      ),
      levels = c("Low", "Intermediate", "High")
  )
  )

data1 <- data0 %>% filter(bs_ab_served_only_as_pre_car_t_bridging_0_no_1_yes==0) ## Exclude patients who received bsAb as bridging

data2 <- data1 %>% 
  filter(best_response %in% c("PR","VGPR","sCR", "CR")) %>%
  filter(!is.na(time_to_best_response)) %>%
  mutate(
    Days.to.death.or.DLC = as.numeric(date_of_last_contact - date_of_bs_ab_step_up_d1) - time_to_best_response,
    Days.to.death.or.relapse.or.DLC = as.numeric(date_progression_or_death  - date_of_bs_ab_step_up_d1)- time_to_best_response,
  ) %>%
  filter(Days.to.death.or.DLC >0 & Days.to.death.or.relapse.or.DLC >0)


## Response
uv_tab <- tbl_uvregression(
  data0 %>%
    dplyr::select(all_of(preds0)),
  method = glm,
  y = `Response`,
  method.args = list(family = binomial, weights = data0$weights),
  exponentiate = TRUE,
  add_estimate_to_reference_rows = FALSE
) %>% modify_header(estimate = "**Estimate**")
## There was a warning constructing the model for variable "BsAb". See message
## below.
## ! non-integer #successes in a binomial glm!
## There was a warning constructing the model for variable "Baseline Hgb". See
## message below.
## ! non-integer #successes in a binomial glm!
## There was a warning constructing the model for variable "Baseline LDH". See
## message below.
## ! non-integer #successes in a binomial glm!
## There was a warning constructing the model for variable "ECOG 2+". See message
## below.
## ! non-integer #successes in a binomial glm!
## There was a warning constructing the model for variable "True EMD". See message
## below.
## ! non-integer #successes in a binomial glm!
## There was a warning constructing the model for variable "Prior LOTs". See
## message below.
## ! non-integer #successes in a binomial glm!
## There was a warning constructing the model for variable "Prior BCMA". See
## message below.
## ! non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "BsAb". See message
## below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "Baseline Hgb". See
## message below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "Baseline LDH". See
## message below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "ECOG 2+". See
## message below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "True EMD". See
## message below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "Prior LOTs". See
## message below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "Prior BCMA". See
## message below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!
mv_tab<-glm( `Response` ~ BsAb + `True EMD` + `ECOG 2+` + `Prior LOTs` + `Prior BCMA` + `Baseline Hgb` + `Baseline LDH`, data0 %>%
    dplyr::select(all_of(preds0), weights = weights) , family = binomial ) %>%
  tbl_regression(exponentiate=TRUE) %>% modify_header(estimate = "**Estimate**") 

table_response <- tbl_merge(
  list(uv_tab, mv_tab),
  tab_spanner = c("**Univariate**", "**Multivariable**")
) 

## VGPR or better
uv_tab <- tbl_uvregression(
  data0 %>%
    dplyr::select(all_of(preds6)),
  method = glm,
  y = `VGPR`,
  method.args = list(family = binomial, weights = data0$weights),
  exponentiate = TRUE,
  add_estimate_to_reference_rows = FALSE,
) %>% modify_header(estimate = "**Estimate**")
## There was a warning constructing the model for variable "BsAb". See message
## below.
## ! non-integer #successes in a binomial glm!
## There was a warning constructing the model for variable "Baseline Hgb". See
## message below.
## ! non-integer #successes in a binomial glm!
## There was a warning constructing the model for variable "Baseline LDH". See
## message below.
## ! non-integer #successes in a binomial glm!
## There was a warning constructing the model for variable "ECOG 2+". See message
## below.
## ! non-integer #successes in a binomial glm!
## There was a warning constructing the model for variable "True EMD". See message
## below.
## ! non-integer #successes in a binomial glm!
## There was a warning constructing the model for variable "Prior LOTs". See
## message below.
## ! non-integer #successes in a binomial glm!
## There was a warning constructing the model for variable "Prior BCMA". See
## message below.
## ! non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "BsAb". See message
## below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "Baseline Hgb". See
## message below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "Baseline LDH". See
## message below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "ECOG 2+". See
## message below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "True EMD". See
## message below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "Prior LOTs". See
## message below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "Prior BCMA". See
## message below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!
mv_tab<-glm( `VGPR` ~ BsAb +`True EMD` + `ECOG 2+` + `Prior LOTs` + `Prior BCMA` + `Baseline Hgb` + `Baseline LDH`, data0 %>%
    dplyr::select(all_of(preds6), weights = weights) , family = binomial ) %>%
  tbl_regression(exponentiate=TRUE) %>% modify_header(estimate = "**Estimate**")

table_VGPR <- tbl_merge(
  list(uv_tab, mv_tab),
  tab_spanner = c("**Univariate**", "**Multivariable**")
) 

## sCR/CR
uv_tab <- tbl_uvregression(
  data0 %>%
    dplyr::select(all_of(preds1)),
  method = glm,
  y = `CR/sCR`,
  method.args = list(family = binomial, weights = data0$weights),
  exponentiate = TRUE,
  add_estimate_to_reference_rows = FALSE,
) %>% modify_header(estimate = "**Estimate**")
## There was a warning constructing the model for variable "BsAb". See message
## below.
## ! non-integer #successes in a binomial glm!
## There was a warning constructing the model for variable "Baseline Hgb". See
## message below.
## ! non-integer #successes in a binomial glm!
## There was a warning constructing the model for variable "Baseline LDH". See
## message below.
## ! non-integer #successes in a binomial glm!
## There was a warning constructing the model for variable "ECOG 2+". See message
## below.
## ! non-integer #successes in a binomial glm!
## There was a warning constructing the model for variable "True EMD". See message
## below.
## ! non-integer #successes in a binomial glm!
## There was a warning constructing the model for variable "Prior LOTs". See
## message below.
## ! non-integer #successes in a binomial glm!
## There was a warning constructing the model for variable "Prior BCMA". See
## message below.
## ! non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "BsAb". See message
## below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "Baseline Hgb". See
## message below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "Baseline LDH". See
## message below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "ECOG 2+". See
## message below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "True EMD". See
## message below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "Prior LOTs". See
## message below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "Prior BCMA". See
## message below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!
mv_tab<-glm( `CR/sCR` ~ BsAb + `True EMD` + `ECOG 2+` + `Prior LOTs` + `Prior BCMA` + `Baseline Hgb` + `Baseline LDH`, data0 %>%
    dplyr::select(all_of(preds1), weights = weights) , family = binomial ) %>%
  tbl_regression(exponentiate=TRUE) %>% modify_header(estimate = "**Estimate**")

table_sCR.CR <- tbl_merge(
  list(uv_tab, mv_tab),
  tab_spanner = c("**Univariate**", "**Multivariable**")
) 

## OS

uv_tab <- tbl_uvregression(
  data1[c(preds2)],
  method = coxph,
  y = Surv(data1$Days.to.death.or.DLC, data1$Death),
  exponentiate = TRUE,
  method.args = list(weights = data1$weights)
) %>% modify_header(estimate = "**Estimate**")

mv_tab<-coxph( Surv(Days.to.death.or.DLC, Death) ~ BsAb + `True EMD` + `ECOG 2+` + `Prior LOTs` + `Prior BCMA` + `Baseline Hgb` + `Baseline LDH`, data1,  weights = weights) %>%
  tbl_regression(exponentiate=TRUE) %>% modify_header(estimate = "**Estimate**")

table_os <- tbl_merge(
  list(uv_tab, mv_tab),
  tab_spanner = c("**Univariate**", "**Multivariable**")
) 


## PFS
uv_tab <- tbl_uvregression(
  data1[c(preds2)],
  method = coxph,
  y = Surv(data1$Days.to.death.or.relapse.or.DLC, data1$Relapse.or.death),
  exponentiate = TRUE,
  method.args = list(weights = data1$weights)
) %>% modify_header(estimate = "**Estimate**")

mv_tab<-coxph( Surv(Days.to.death.or.relapse.or.DLC, Relapse.or.death) ~ BsAb + `True EMD` + `ECOG 2+` + `Prior LOTs` + `Prior BCMA` + `Baseline Hgb` + `Baseline LDH`, weights = weights, data1) %>%
  tbl_regression(exponentiate=TRUE) %>% modify_header(estimate = "**Estimate**")

table_pfs <- tbl_merge(
  list(uv_tab, mv_tab),
  tab_spanner = c("**Univariate**", "**Multivariable**")
)

## DOR
uv_tab <- tbl_uvregression(
  data2[c(preds2)],
  method = coxph,
  y = Surv(data2$Days.to.death.or.relapse.or.DLC, data2$Relapse.or.death),
  exponentiate = TRUE,
  method.args = list(weights = data2$weights)
) %>% modify_header(estimate = "**Estimate**")

mv_tab<-coxph( Surv(  Days.to.death.or.relapse.or.DLC, Relapse.or.death) ~ BsAb + `True EMD` + `ECOG 2+` + `Prior LOTs` + `Prior BCMA` + `Baseline Hgb` + `Baseline LDH`, weights = weights, data2) %>%
  tbl_regression(exponentiate=TRUE) %>% modify_header(estimate = "**Estimate**")

table_dor <- tbl_merge(
  list(uv_tab, mv_tab),
  tab_spanner = c("**Univariate**", "**Multivariable**")
)

## IFS
uv_tab <- tbl_uvregression(
  data1[c(preds2)],
  method = coxph,
  y = Surv(data1$days_to_infection_or_last_contact, data1$infection_or_death),
  exponentiate = TRUE,
  method.args = list(weights = data1$weights)
 ) %>% modify_header(estimate = "**Estimate**")

mv_tab<-coxph( Surv(  days_to_infection_or_last_contact, infection_or_death) ~ BsAb + `True EMD` + `ECOG 2+` + `Prior LOTs` + `Prior BCMA` + `Baseline Hgb` + `Baseline LDH`, weights = weights, data1) %>%
  tbl_regression(exponentiate=TRUE) %>% modify_header(estimate = "**Estimate**")

table_ifs <- tbl_merge(
  list(uv_tab, mv_tab),
  tab_spanner = c("**Univariate**", "**Multivariable**")
)

## ICANS 1+
uv_tab <- tbl_uvregression(
  data0 %>%
    dplyr::select(all_of(preds3)),
  method = glm,
  y = `ICANS 1+`,
  method.args = list(family = binomial, weights = data0$weights),
  exponentiate = TRUE,
  add_estimate_to_reference_rows = FALSE
) %>% modify_header(estimate = "**Estimate**")
## There was a warning constructing the model for variable "BsAb". See message
## below.
## ! non-integer #successes in a binomial glm!
## There was a warning constructing the model for variable "Baseline Hgb". See
## message below.
## ! non-integer #successes in a binomial glm!
## There was a warning constructing the model for variable "Baseline LDH". See
## message below.
## ! non-integer #successes in a binomial glm!
## There was a warning constructing the model for variable "ECOG 2+". See message
## below.
## ! non-integer #successes in a binomial glm!
## There was a warning constructing the model for variable "True EMD". See message
## below.
## ! non-integer #successes in a binomial glm!
## There was a warning constructing the model for variable "Prior LOTs". See
## message below.
## ! non-integer #successes in a binomial glm!
## There was a warning constructing the model for variable "Prior BCMA". See
## message below.
## ! non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "BsAb". See message
## below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "Baseline Hgb". See
## message below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "Baseline LDH". See
## message below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "ECOG 2+". See
## message below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "True EMD". See
## message below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "Prior LOTs". See
## message below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "Prior BCMA". See
## message below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!
mv_tab<-glm( `ICANS 1+` ~ BsAb + `True EMD` + `ECOG 2+` + `Prior LOTs` + `Prior BCMA` + `Baseline Hgb` + `Baseline LDH`, data0 %>%
    dplyr::select(all_of(preds3), weights = weights) , family = binomial ) %>%
  tbl_regression(exponentiate=TRUE) %>% modify_header(estimate = "**Estimate**")

table_ICANS <- tbl_merge(
  list(uv_tab, mv_tab),
  tab_spanner = c("**Univariate**", "**Multivariable**")
) 

## CRS 1+
uv_tab <- tbl_uvregression(
  data0 %>%
    dplyr::select(all_of(preds4)),
  method = glm,
  y = `CRS 1+`,
  method.args = list(family = binomial, weights = data0$weights),
  exponentiate = TRUE,
  add_estimate_to_reference_rows = FALSE
) %>% modify_header(estimate = "**Estimate**")
## There was a warning constructing the model for variable "BsAb". See message
## below.
## ! non-integer #successes in a binomial glm!
## There was a warning constructing the model for variable "Baseline Hgb". See
## message below.
## ! non-integer #successes in a binomial glm!
## There was a warning constructing the model for variable "Baseline LDH". See
## message below.
## ! non-integer #successes in a binomial glm!
## There was a warning constructing the model for variable "ECOG 2+". See message
## below.
## ! non-integer #successes in a binomial glm!
## There was a warning constructing the model for variable "True EMD". See message
## below.
## ! non-integer #successes in a binomial glm!
## There was a warning constructing the model for variable "Prior LOTs". See
## message below.
## ! non-integer #successes in a binomial glm!
## There was a warning constructing the model for variable "Prior BCMA". See
## message below.
## ! non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "BsAb". See message
## below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "Baseline Hgb". See
## message below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "Baseline LDH". See
## message below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "ECOG 2+". See
## message below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "True EMD". See
## message below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "Prior LOTs". See
## message below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "Prior BCMA". See
## message below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!
mv_tab<-glm( `CRS 1+` ~ BsAb + `True EMD` + `ECOG 2+` + `Prior LOTs` + `Prior BCMA` + `Baseline Hgb` + `Baseline LDH`, data0 %>%
    dplyr::select(all_of(preds4), weights = weights) , family = binomial ) %>%
  tbl_regression(exponentiate=TRUE) %>% modify_header(estimate = "**Estimate**")

table_CRS <- tbl_merge(
  list(uv_tab, mv_tab),
  tab_spanner = c("**Univariate**", "**Multivariable**")
) 

## Steroid use
uv_tab <- tbl_uvregression(
  data0 %>%
    dplyr::select(all_of(preds5)),
  method = glm,
  y = `Steroid use`,
  method.args = list(family = binomial, weights = data0$weights),
  exponentiate = TRUE,
  add_estimate_to_reference_rows = FALSE
) %>% modify_header(estimate = "**Estimate**")
## There was a warning constructing the model for variable "BsAb". See message
## below.
## ! non-integer #successes in a binomial glm!
## There was a warning constructing the model for variable "Baseline Hgb". See
## message below.
## ! non-integer #successes in a binomial glm!
## There was a warning constructing the model for variable "Baseline LDH". See
## message below.
## ! non-integer #successes in a binomial glm!
## There was a warning constructing the model for variable "ECOG 2+". See message
## below.
## ! non-integer #successes in a binomial glm!
## There was a warning constructing the model for variable "True EMD". See message
## below.
## ! non-integer #successes in a binomial glm!
## There was a warning constructing the model for variable "Prior LOTs". See
## message below.
## ! non-integer #successes in a binomial glm!
## There was a warning constructing the model for variable "Prior BCMA". See
## message below.
## ! non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "BsAb". See message
## below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "Baseline Hgb". See
## message below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "Baseline LDH". See
## message below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "ECOG 2+". See
## message below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "True EMD". See
## message below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "Prior LOTs". See
## message below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!
## There was a warning running `tbl_regression()` for variable "Prior BCMA". See
## message below.
## ! non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, non-integer
##   #successes in a binomial glm!, non-integer #successes in a binomial glm!,
##   non-integer #successes in a binomial glm!, non-integer #successes in a
##   binomial glm!, non-integer #successes in a binomial glm!, …, non-integer
##   #successes in a binomial glm!, and non-integer #successes in a binomial glm!
mv_tab<-glm( `Steroid use` ~ BsAb + `True EMD` + `ECOG 2+` + `Prior LOTs` + `Prior BCMA` + `Baseline Hgb` + `Baseline LDH`, data0 %>%
    dplyr::select(all_of(preds5), weights = weights) , family = binomial ) %>%
  tbl_regression(exponentiate=TRUE) %>% modify_header(estimate = "**Estimate**")

table_steroid <- tbl_merge(
  list(uv_tab, mv_tab),
  tab_spanner = c("**Univariate**", "**Multivariable**")
) 


# Combine all tables into one
final_table <- tbl_stack(
  list(
    table_response %>% modify_header() %>% modify_footnote(everything() ~ NA, abbreviation = TRUE) %>%
      modify_footnote(estimate_1 ~ "Estimate = Odds ratio or hazard ratio", abbreviation = TRUE),
    table_VGPR,
    table_sCR.CR,
    table_os,
    table_pfs,
    table_dor,
    table_ifs,
    table_ICANS,
    table_CRS,
    table_steroid
  ),
  group_header = c(
    "Response",
    "≥VGPR",
    "CR/sCR",
    "Overall Survival",
    "Progression-free Survival",
    "Duration of Response",
    "Infection-free Survival",
    "ICANS 1+",
    "CRS 1+",
    "Steroid use"
  )
)
  
final_table
Characteristic
Univariate
Multivariable
N Estimate1 95% CI p-value Estimate 95% CI p-value
Response
BsAb 370





    Tec
— —
— —
    Elra
0.86 0.55, 1.35 0.5 0.86 0.52, 1.42 0.6
Baseline Hgb 370 1.46 1.29, 1.65 <0.001 1.34 1.18, 1.53 <0.001
Baseline LDH 370 0.86 0.77, 0.94 0.002 0.86 0.77, 0.95 0.004
ECOG 2+ 370 0.39 0.25, 0.61 <0.001 0.42 0.26, 0.69 <0.001
True EMD 370 0.91 0.52, 1.63 0.8 0.94 0.50, 1.78 0.9
Prior LOTs 370 0.87 0.79, 0.94 0.001 0.91 0.82, 1.01 0.087
Prior BCMA 370 0.56 0.36, 0.85 0.008 0.50 0.29, 0.87 0.014
≥VGPR
BsAb 370





    Tec
— —
— —
    Elra
0.62 0.40, 0.95 0.031 0.57 0.35, 0.93 0.026
Baseline Hgb 370 1.38 1.23, 1.55 <0.001 1.32 1.17, 1.49 <0.001
Baseline LDH 370 0.84 0.73, 0.93 0.003 0.84 0.73, 0.93 0.004
ECOG 2+ 370 0.52 0.33, 0.81 0.004 0.57 0.34, 0.93 0.026
True EMD 370 0.63 0.35, 1.10 0.11 0.56 0.30, 1.04 0.071
Prior LOTs 370 0.88 0.80, 0.95 0.003 0.95 0.86, 1.06 0.4
Prior BCMA 370 0.47 0.31, 0.72 <0.001 0.41 0.24, 0.70 0.001
CR/sCR
BsAb 370





    Tec
— —
— —
    Elra
0.70 0.44, 1.10 0.12 0.68 0.40, 1.12 0.13
Baseline Hgb 370 1.44 1.29, 1.63 <0.001 1.39 1.23, 1.58 <0.001
Baseline LDH 370 0.86 0.74, 0.96 0.019 0.87 0.75, 0.98 0.044
ECOG 2+ 370 0.51 0.31, 0.81 0.006 0.61 0.36, 1.03 0.068
True EMD 370 0.66 0.35, 1.20 0.2 0.62 0.31, 1.19 0.2
Prior LOTs 370 0.89 0.81, 0.97 0.014 0.95 0.85, 1.06 0.4
Prior BCMA 370 0.56 0.37, 0.86 0.009 0.52 0.31, 0.89 0.018
Overall Survival
BsAb 352





    Tec
— —
— —
    Elra
1.41 0.96, 2.07 0.083 1.46 0.97, 2.20 0.073
Baseline Hgb 352 0.73 0.66, 0.81 <0.001 0.77 0.69, 0.86 <0.001
Baseline LDH 352 1.13 1.06, 1.20 <0.001 1.15 1.10, 1.20 <0.001
ECOG 2+ 352 2.25 1.61, 3.16 <0.001 1.92 1.32, 2.78 <0.001
True EMD 352 1.34 0.88, 2.03 0.2 1.28 0.82, 2.00 0.3
Prior LOTs 352 1.06 0.99, 1.13 0.091 1.04 0.97, 1.11 0.3
Prior BCMA 352 1.09 0.77, 1.54 0.6 1.20 0.81, 1.78 0.4
Progression-free Survival
BsAb 355





    Tec
— —
— —
    Elra
1.43 1.05, 1.94 0.021 1.52 1.11, 2.09 0.009
Baseline Hgb 355 0.77 0.71, 0.83 <0.001 0.81 0.74, 0.89 <0.001
Baseline LDH 355 1.18 1.08, 1.28 <0.001 1.17 1.08, 1.27 <0.001
ECOG 2+ 355 1.63 1.23, 2.16 <0.001 1.43 1.04, 1.95 0.027
True EMD 355 1.04 0.74, 1.47 0.8 0.94 0.67, 1.32 0.7
Prior LOTs 355 1.10 1.05, 1.16 <0.001 1.06 1.00, 1.12 0.056
Prior BCMA 355 1.58 1.18, 2.10 0.002 1.63 1.17, 2.27 0.004
Duration of Response
BsAb 177





    Tec
— —
— —
    Elra
1.49 0.87, 2.56 0.14 1.50 0.85, 2.63 0.2
Baseline Hgb 177 0.84 0.74, 0.95 0.007 0.83 0.72, 0.96 0.014
Baseline LDH 177 1.05 0.90, 1.23 0.5 1.04 0.89, 1.21 0.7
ECOG 2+ 177 1.05 0.64, 1.72 0.8 0.89 0.52, 1.52 0.7
True EMD 177 1.48 0.75, 2.93 0.3 1.36 0.68, 2.72 0.4
Prior LOTs 177 1.01 0.89, 1.15 0.8 1.02 0.90, 1.17 0.7
Prior BCMA 177 1.08 0.67, 1.74 0.8 1.06 0.65, 1.75 0.8
Infection-free Survival
BsAb 358





    Tec
— —
— —
    Elra
1.25 0.93, 1.67 0.13 1.24 0.92, 1.67 0.2
Baseline Hgb 358 0.87 0.81, 0.94 <0.001 0.90 0.83, 0.97 0.010
Baseline LDH 358 1.10 1.05, 1.15 <0.001 1.09 1.05, 1.14 <0.001
ECOG 2+ 358 1.70 1.30, 2.22 <0.001 1.52 1.14, 2.02 0.004
True EMD 358 1.16 0.83, 1.62 0.4 1.10 0.79, 1.54 0.6
Prior LOTs 358 0.98 0.93, 1.03 0.4 0.98 0.93, 1.04 0.5
Prior BCMA 358 0.85 0.65, 1.11 0.2 0.99 0.73, 1.35 >0.9
ICANS 1+
BsAb 370





    Tec
— —
— —
    Elra
0.87 0.49, 1.52 0.6 0.91 0.49, 1.65 0.8
Baseline Hgb 370 0.83 0.72, 0.95 0.008 0.84 0.72, 0.98 0.024
Baseline LDH 370 1.05 0.95, 1.15 0.3 1.04 0.93, 1.14 0.5
ECOG 2+ 370 3.03 1.77, 5.21 <0.001 2.44 1.39, 4.31 0.002
True EMD 370 1.40 0.69, 2.69 0.3 1.31 0.62, 2.62 0.5
Prior LOTs 370 0.88 0.77, 0.98 0.034 0.88 0.76, 1.01 0.069
Prior BCMA 370 0.61 0.36, 1.03 0.067 1.01 0.53, 1.92 >0.9
CRS 1+
BsAb 370





    Tec
— —
— —
    Elra
0.53 0.34, 0.82 0.004 0.55 0.35, 0.86 0.010
Baseline Hgb 370 0.94 0.85, 1.04 0.3 0.91 0.81, 1.01 0.076
Baseline LDH 370 0.98 0.90, 1.06 0.5 0.97 0.89, 1.06 0.5
ECOG 2+ 370 1.08 0.70, 1.68 0.7 0.97 0.60, 1.54 0.9
True EMD 370 0.58 0.33, 1.02 0.061 0.53 0.29, 0.94 0.033
Prior LOTs 370 0.92 0.84, 1.00 0.045 0.87 0.79, 0.96 0.006
Prior BCMA 370 1.17 0.77, 1.76 0.5 1.84 1.12, 3.04 0.016
Steroid use
BsAb 355





    Tec
— —
— —
    Elra
0.83 0.48, 1.40 0.5 0.87 0.50, 1.50 0.6
Baseline Hgb 355 0.89 0.79, 1.01 0.074 0.90 0.79, 1.03 0.14
Baseline LDH 355 1.07 0.98, 1.17 0.12 1.05 0.96, 1.15 0.3
ECOG 2+ 355 1.71 1.03, 2.83 0.036 1.32 0.78, 2.23 0.3
True EMD 355 1.61 0.85, 2.97 0.13 1.57 0.81, 2.94 0.2
Prior LOTs 355 0.88 0.79, 0.98 0.028 0.94 0.83, 1.06 0.3
Prior BCMA 355 0.47 0.29, 0.78 0.003 0.65 0.36, 1.16 0.15
1 Estimate = Odds ratio or hazard ratio

8 CONTOUR PLOTS

8.1 FIGURE 4: Baseline Hgb and LDH

data1 <- IPTW.data %>%
  filter(bs_ab_served_only_as_pre_car_t_bridging_0_no_1_yes==0) %>% ## Exclude patients who received bsAb as bridging
  #filter(bs_ab_name == "tec") %>%
  transmute(
    bs_ab_name,
    weights,
    "CR/sCR" = ifelse(best_response %in% c("CR", "sCR"),1,0),
    "PD" = ifelse(best_response %in% c("PD"),1,0),
    "True EMD" = ifelse(just_prior_to_bs_ab_initiation_is_there_true_emd_if_none_move_to_column_gk == 2, 1,0),
    "Age" = age_at_bs_ab_initiation/10,
    "ECOG" = ecog_ps_at_bs_ab_initiation_0_5,
    ECOG_2_4 = case_when(
      is.na(ecog_ps_at_bs_ab_initiation_0_5) ~ NA,
      ecog_ps_at_bs_ab_initiation_0_5 %in% c(2,3,4) ~ 1,
      TRUE ~ 0),
    "Prior LOTs" = number_of_prior_lines_of_therapy_before_bs_ab,
    "HRCA" = ifelse(deletion_17p_prior_to_bs_ab_0_no_1_yes | t_4_14_prior_to_bs_ab_0_no_1_yes | t_14_16_prior_to_bs_ab_0_no_1_yes | gain_amp1q21_prior_to_bs_ab_0_no_1_yes, 1, 0),
    "Prior BCMA" = prior_bcma_directed_therapy_yes_no,
    Relapse.or.death = ifelse(progression_yes_no | death_yes_no,1,0),
    Death = ifelse(death_yes_no,1,0),
    date_of_pd,
    date_of_last_contact,
    date_progression_or_death = coalesce(date_of_pd, date_of_last_contact),
    Days.to.death.or.DLC = as.numeric(date_of_last_contact - date_of_bs_ab_step_up_d1),
    Days.to.death.or.relapse.or.DLC = as.numeric(date_progression_or_death  - date_of_bs_ab_step_up_d1),
    "ALC BL" = baseline_alc,
    "CrCl BL" = cr_cl_prior_to_step_up_dose_number_1,
    "Trial ineligible" = ifelse(eligible_for_corresponding_bs_ab_trial_on_first_dose_use_trial_elig_sheet_to_help_you == 1,0,1),
    "Baseline platelets" = baseline_pl_ts,
    baseline_hgb,
    baseline_ldh = as.double(baseline_ldh_prior_to_bs_ab_initiation),
    baseline_ferritin = as.numeric(baseline_ferritin_prior_to_bs_ab_initiation),
    ) %>% filter(!is.na(baseline_ldh))


## Tec
data_tec = data1 %>% filter(bs_ab_name == "tec")

dd_tec <- datadist(data_tec)
options(datadist='dd_tec')

S_tec <- Surv(data_tec$Days.to.death.or.relapse.or.DLC/30, data_tec$Relapse.or.death)
f_tec <- cph(S_tec ~ rcs(baseline_hgb, 3) * rcs(baseline_ldh, 3) , x=TRUE, y=TRUE,surv=TRUE,data=data_tec, weights=weights)
med_tec <- Quantile(f_tec)

pred_tec  <- data.frame(Predict(f_tec,  baseline_hgb, baseline_ldh,
                                fun = function(x) med_tec(lp = x))) %>%
  mutate(drug = "Teclistamab")

## Elra
data_elra = data1 %>% filter(bs_ab_name == "elra")

dd_elra <- datadist(data_elra)
options(datadist='dd_elra')

S_elra <- Surv(data_elra$Days.to.death.or.relapse.or.DLC/30, data_elra$Relapse.or.death)
f_elra <- cph(S_elra ~ rcs(baseline_hgb, 3) * rcs(baseline_ldh, 3) , x=TRUE, y=TRUE,surv=TRUE,data=data_elra, weights=weights)
med_elra <- Quantile(f_elra)

pred_elra <- data.frame(Predict(f_elra, baseline_hgb, baseline_ldh,
                                fun = function(x) med_elra(lp = x))) %>%
  mutate(drug = "Elranatamab")

## Combine both predictions
all_pred <- bind_rows(pred_tec, pred_elra)
z_range <- range(all_pred$yhat, na.rm = TRUE)
brks    <- pretty(z_range, n = 10)   # tweak n if you want more/fewer bands

all_pred <- all_pred %>%
  mutate(
    z_band = cut(
      yhat,
      breaks        = brks,
      include.lowest = TRUE,
      right          = FALSE
    )
  )

z_levels <- levels(all_pred$z_band)

cols <- viridis_pal()(length(z_levels))


contour_plot_tec <- all_pred %>%
  filter(drug == "Teclistamab") %>%
  ggplot(aes(baseline_hgb, baseline_ldh, fill = z_band)) +
  geom_raster() +   # or geom_tile(); blocky but guarantees alignment
  theme_classic() +
  theme(
    panel.grid.major = element_line(color = "grey85"),
    panel.grid.minor = element_line(color = "grey92")
  )+
  labs(
    title = "A) Teclistamab",
    x = "Baseline Hgb (g/dL)",
    y = "Baseline LDH (U/L)",
    fill = "Median PFS (months)"
  ) +
  scale_fill_manual(
    values       = cols,
    breaks       = z_levels,
    drop         = FALSE,
    na.translate = FALSE
  ) +
  scale_x_continuous(breaks = seq(0, 20, 1))+
  scale_y_continuous(breaks = seq(0, 2000, 100))+
  coord_cartesian(xlim = c(7, 13), ylim = c(120, 780))

contour_plot_tec

contour_plot_elra <- all_pred %>%
  filter(drug == "Elranatamab") %>%
  ggplot(aes(baseline_hgb, baseline_ldh, fill = z_band)) +
  geom_raster() +
  theme_classic() +
  theme(
    legend.position = "none",
    panel.grid.major = element_line(color = "grey85"),
    panel.grid.minor = element_line(color = "grey92")
    ) +
  labs(
    title = "B) Elranatamab",
    x = "Baseline Hgb (g/dL)",
    y = "Baseline LDH (U/L)"
    #fill = "Median PFS (months)"
  ) +
  scale_fill_manual(
    values       = cols,
    breaks       = z_levels,
    drop         = FALSE,
    na.translate = FALSE
  ) +
  scale_x_continuous(breaks = seq(0, 20, 1))+
  scale_y_continuous(breaks = seq(0, 2000, 100))+
  coord_cartesian(xlim = c(7, 13), ylim = c(120, 780))

 contour_plot_elra

8.2 ***FIGURE 4: Hgb and LDH contour

combined_contour <- (contour_plot_tec / contour_plot_elra) +
  plot_layout(guides = "collect", axes = "collect")


ggsave("svg/FIGURE4.svg", combined_contour, device = "svg",
       width = 6.5,
       height = 6.5,
       units = c("in"))
## Warning: Removed 37 rows containing missing values or values outside the scale range
## (`geom_raster()`).
## Warning: Removed 4293 rows containing missing values or values outside the scale range
## (`geom_raster()`).
ggsave("jpeg/FIGURE4.jpeg", combined_contour, device = "jpeg",
       width = 6.5,
       height = 6.5,
       dpi = 600,
       units = c("in"))
## Warning: Removed 37 rows containing missing values or values outside the scale range
## (`geom_raster()`).
## Removed 4293 rows containing missing values or values outside the scale range
## (`geom_raster()`).
ggsave("tiff/FIGURE4.tiff", combined_contour, device = "tiff",
       width = 6.5,
       height = 6.5,
       dpi = 600,
       units = c("in"),
       compression = "lzw")
## Warning: Removed 37 rows containing missing values or values outside the scale range
## (`geom_raster()`).
## Removed 4293 rows containing missing values or values outside the scale range
## (`geom_raster()`).
knitr::include_graphics("jpeg/FIGURE4.jpeg")